問題タブ [runge-kutta]

C で 4 次のルンゲクッタを使用して、惑星軌道の微分方程式を統合しようとしています。最初に統合しようとした方程式は、位置 dr/dt = sqrt((2/mu)*(E-( k/r)) - pow(l, 2)/(pow(mu, 2)*pow(r, 2)))。プログラムは正しくコンパイルされますが、nan を返し続けます。私が使用しているソースコードは次のとおりです。

私のコードは次のとおりです。

x0 = 1.0、t0 = 0.0 の Rk4 関数を dt の関数としてプロットしようとしています。関数をベクトル化し、タイムステップ dt の配列を作成しようとしましたが、「ValueError: シーケンスで配列要素を設定しています」というエラーが発生します。

適応ステップ サイズで 4 次ルンゲ クッタ数値近似を実行するコードを書いています。

ただし、コードを実行すると、次のように出力されます。

したがって、エラー E はどちらの場合も 10^-5 未満ですが、x(1) は出力されません。

Runge-Kutta4 を使用して、ODE ソリューション用に次のコードを作成しました。

fは関数y' = f(x,y)、y0は初期値、x0は関数の開始位置、hサブインターバル、x_finalは関数の停止位置です。

コードを試してみたところ、ODE が正しく解決されましたが、サブインターバルx0を使用して、間隔の xy 軸上に関数をプロットしたいと考えています。を使用してプロットしようとすると、空のグラフしか表示されません。プロットするために複数にバインドする必要があることを（推測して）理解していますが、コードをあまり変更せずにそれを行うにはどうすればよいですか?x_finalhplot(x0:h:x_final,y)yx

ygiven y0、 interval x0to x_final、 givenのグラフをプロットするにはどうすればよいhですか?

MATLAB は初めてなので、できる限りのサポートをお願いします!

編集:私のコードの目的を明確にするため。

ソリューションとグラフ作成の両方に、この ODE ソルバーが必要です。yonの値を とh比較して切り捨て誤差を研究し、 のグラフを異なるで比較し2*hて Runge-Kutta4 の安定性を研究することになっています。 yh

結合された DE を使用して、Lotka-Volterra Predator-Prey システムをモデル化しようとしています。

dy(1)/dt = rx(1-x/k) - ay(1)y(2) % 獲物の個体数

dy(2)/dt = aby(1)y(2) - dy(2) % 捕食者個体群

これが私が持っているコードです:

時系列プロット (y(1) 対 時間) が定常状態ではなく振動を示すようにコードを変更するにはどうすればよいですか? プロットは私が望むように構成されていますが、関数の動作は私が期待したものではないため、統合ステップで何かがうまくいかないようです。

定常状態の時系列プロット

平衡に達するフェーズ プロット