問題タブ [runge-kutta]

この微分方程式を解くように求められました。

(x,y,vx,vy)'=(vx,vy,vy,-vx)

2*piピリオド付きの円運動を返す必要があります。私は機能を実装しました：

に_alpha等しい場所0。さて、これはオイラーの方法で問題なく機能します。この ODE を に対して積分すると2 * pi * 10、初期条件 が与えられ、(1, 0, 0, -1)ある精度で、予想どおり、位置がの範囲内で に0.003匹敵することがわかります。しかし、同じ ODE を次のように実装されたルンゲ・クッタの方法 (精度で数秒) と統合すると、次のようになります。(1, 0)1 ± 0.10.0032 * pi * 10

プログラムは、4 次のルンゲ クッタ法で、理論的には精度が必要な場合に、およそ にx等しい位置を返します。確認したところ、Runge_Kutta の周期はほぼ 4 倍になっているように見えます (経過中に からに到達するため) が、その理由がわかりません。これは私のメインの内容です：0.391E-62 * pix10.48

前もって感謝します。

Lotka-Volterra systtem を解決するために runge-kutta メソッドを実装しようとしていますが、コード (以下) が正しく機能していません。StackOverflow の他のトピックで見つけた推奨事項に従いましたが、結果は、たとえば Pylab で利用可能な rk4 メソッドのような組み込みの Runge-Kutta メソッドでは収束しません。誰かが私を助けることができますか？

コメントに従ってコードを更新しました。したがって、関数meurk4：

になりました（修正済み）：

それにもかかわらず、結果は次のとおりです。

ここに画像の説明を入力

組み込みメソッド rk4 (Pylab から) の結果は次のようになります。

ここに画像の説明を入力

したがって、結果が組み込みの rk4 メソッドと同じではないため、確かに私のコードはまだ正しくありません。誰か助けてください。

このコードは大丈夫ですか？

以下に示す時間テストでは、Skyfieldobj.at(jd).position.kmがの単一の時間値を返すのに数百マイクロ秒から 1 ミリ秒かかることがわかりましたjdが、より長いオブジェクト (時間内のポイントのリスト) の増分コストはJulianDate、ポイントごとに約 1 マイクロ秒に過ぎません。 . Jplephemと 2 つの異なるエフェメリスを使用した場合、同様の速度が見られます。

ここでの私の質問は次のとおりです。たとえば、独自の可変ステップサイズを使用する外部ルンゲクッタルーチンのスレーブとして、特定の時点にランダムアクセスしたい場合、Python内でこれをより高速に実行できる方法はありますか(コードのコンパイル方法を学ぶ)?

これは、Skyfield の一般的な使用方法ではないことを理解しています。通常、JulianDateオブジェクトを長いリストの時点でロードし、それらを一度に計算します。軌道積分器が行う方法のように、おそらく数千回 (またはそれ以上) ではなく、数回実行します。

回避策:時間粒度の細かいオブジェクトをNumPy使用して Skyfield を 1 回実行し、タイムステップが常にJulianDateNumPy 配列のストライディングに直接対応します。

または、再補間を試すこともできます。私は非常に正確な計算を行っていないので、単純な NumPy または SciPy の 2 次で問題ないかもしれません。

最終的には、太陽系の重力場の影響下にあるオブジェクト (深宇宙衛星、彗星、小惑星など) の経路を統合してみたいと思います。軌道解を探すとき、6D 位相空間で何百万もの開始状態ベクトルを試すかもしれません。ob.at(jd).observe(large_body).position.km重力は他のすべてのものと同じように光速で移動するため、メソッドなどを使用する必要があることはわかっています。これは反復計算であるため（推測では）かなりの時間がかかるようです（「見てみましょう...木星はどこにあったので、今すぐ重力だと感じます」）。しかし、コズミックオニオンを一度に1層ずつ剥がしましょう.

図 1. de405 と de421 のさまざまな長さJulianDateのオブジェクトに対する私のラップトップでの Skyfield と JPLephem のパフォーマンス。それらはすべてほぼ同じです-（非常に）最初のポイントで約0.5ミリ秒、追加のポイントごとに1マイクロ秒です。また、スクリプトの実行時に計算される最初のポイント(Earth (blue) with len(jd) = 1) には、追加のミリ秒のアーティファクトがあります。

地球と月は、内部で 2 段階の計算 (地球と月の重心と重心に関する個々の軌道) であるため、処理が遅くなります。水星は、エフェメリスの時間ステップと比較して非常に速く移動するため、（コストのかかる）チェビシェフ補間でより多くの係数が必要になるため、遅くなる可能性がありますか?

SCRIPT FOR SKYFIELD DATA JPLephem スクリプトはさらに下にあります

SCRIPT FOR JPLEPHEM DATA Skyfield スクリプトは上にあります

fdm_2nd とガウス肉屋係数を使用して、RK 暗黙的な 2 次から対流拡散方程式 (1D) を実装しようとしています: 'u_t = -uu_x + nu .u_xx' 。

私の目標は、明示的スキームと暗黙的スキームを比較することです。粘度の少ない数値でうまく機能する明示的なrk。明示的なスキームの曲線は、非常に素晴らしい衝撃波を示しています。

k(i) 係数のソルバーを正しく実装するには、あなたの助けが必要です。すべての k(i) に対して newton メソッドを実装する方法がわかりません。すべての時空間ステップに実装する必要がありますか? またはジャストインタイム？ヤコビアンは間違っているかもしれませんが、どこにあるのかわかりません。それとも、ジャコビアンを間違った方向に使用しているのかもしれません...

実際、私のコードは機能しますが、どこか間違っていると思います...また、暗黙の曲線は初期値から移動しません。

ここに私の機能：

ヤコビアン :

ここに私のニュートンコード：

4次のルンゲ・クッタ積分を使用して振り子の速度と力を計算するコードを書きましたが、残念ながら次のエラーが発生するため実行できません。

何が欠けているのか、何が間違いなのかわかりません。どんなアイデアでも役に立ちます。ありがとう！

これは私のコードです:

私は、4 つのリンクされた微分方程式のオイラー法実装を 4 次ルンゲ クッタ実装に変換してきました。私は一般的なアプローチが正しいと合理的に確信しており、RK4の適用方法を理解していますが、おそらく6年間、中途半端な深刻な数学を行っていないため、何かが欠けている可能性があります. 私の RK4 計算は、ステップ サイズ 1 を使用すると適切な出力が得られますが、他のステップ サイズを使用するとすぐにゼロになります。新鮮な目が私が間違ったことをすぐに見つけられることを願っています. 完全な解決策を投稿しないでください-私が作った可能性のあるエラーへのポインタを好むでしょう-それがコードであろうとRK4の理解であろうと、私はこれを自分で理解できるようにしたいからです。

これが私の Euler の実装です。うまく機能します

これは私のRK4実装です

これは速度の問題です。次の 3 つの動作状態を持つ拡散方程式を解こうとしています。

• ラムダ == 0 平衡
• ラムダ > 0 最大拡散
• ラムダ < 0 分拡散

ボトルネックは、拡散演算子関数のelseステートメントです

平衡状態には、単純な T 演算子と拡散演算子があります。他の 2 つの州ではかなり複雑になります。これまでのところ、コードの実行時間をじっと待つ忍耐力がありませんでした。方程式が正しく、平衡状態の出力が正しいように見える限り、おそらく誰かが非平衡状態の速度を上げるためのヒントを持っていますか?

(Runge-Kutta の代わりに Euler 解 (FTCS) の方が速いと思います。これはまだ試していません。)

白黒の画像をインポートして、コードを試すことができます。

要するに、f_out_equil は計算が比較的高速ですが、min および max ケースは高価で時間がかかります。

私が使用している画像へのリンクは次のとおりです。

コーディングを改善するためのヒントを感謝します, どうもありがとう,

出力用の簡単なプロット スクリプトを次に示します。