問題タブ [runge-kutta]

Verlet 積分が Euler 積分より優れている理由を誰か説明してもらえますか? そして、なぜ RK4 は Verlet よりも優れているのでしょうか? なぜそれがより良い方法なのか理解できません。

この引用は主にこのスレッドの結果です：Javaの微分方程式。
基本的に、私はJason S.のアドバイスに従い、ルンゲクッタ法（RK4）を介して微分方程式の数値解法を実装しようとしました。

みなさん、こんにちは。JavaでSIRエピデミックモデルの簡単なシミュレーションプログラムを作成しようとしています。基本的に、SIRは次の3つの微分方程式のシステムによって定義されます
。S'（t）= --lamda（t）* S（t）
I'（t）= lamda（t）* S（t）-gamma（t）* I（t）
R'（t）= gamma（t）* I（t）
S-影響を受けやすい人々、I-感染した人々、R-回復した人々。lamda（t）= [c * x * I（t）] / N（T）c-接触の数、x-感染性（病気の人との接触後に病気になる確率）、N（t）-総人口（一定です）。
gamma（t）= 1 /病気の期間（一定）

最初の試みはあまり成功しませんでしたが、ルンゲクッタ法でこの方程式を解こうとしましたが、この試みの結果、次のコードが生成されました。

病気の人（I）の数は最初に増加し、次に約0に減少し、回復した人の数は厳密に増加する必要があるため、これは機能していないようです。病気+健康+回復の総数は100になるはずですが、私のコードはいくつかの奇妙な結果を生成します：

間違いは見つかりません、アドバイスしてください！よろしくお願いします！

Runge-Kuttaの4次積分法を2次元行列（x、y）のmatlabに実装しようとしています。マトリックスには、各ポイントの高さの値（float）が含まれています。アイデアは、マトリックスにパーティクルを配置し、方向ベクトルuおよびvによって「プッシュ」されるときのパスを監視することです。

matlabに組み込まれている「gradient」関数を使用して、オイラーの積分手法を実装しました。ただし、RK4統合の場合、勾配によって与えられる所定の勾配ではなく、可変点で4つの勾配を取る必要があります。

Matlabの勾配関数自体にはこの種のパラメーターがないと思います（Matlabのドキュメントから理解していることから）。これにはどのようなアプローチを取るべきですか？

任意の点で勾配を計算する独自の勾配関数を作成する場合、それは直進関数ではないため、どのように作成するかはよくわかりません。このアプローチについて何か提案はありますか？

ありがとう。

ルンゲクッタ法（「RK4」）が常微分方程式の正確な解の0.01％以内になるために必要なステップ数を決定しようとしています。これをオイラー法と比較しています。どちらも、両対数プロット上で直線になるはずです。私のオイラーソリューションは正しいようですが、RKのカーブしたソリューションを取得しています。それらは同じコードに基づいているので、私は問題について完全に混乱しています。

編集：ウィキペディアのリンクを削除して申し訳ありません。私は新しいユーザーなので、複数のリンクを保持することはできません。ただし、どちらの方法も、私の実装と同様にWikipediaで詳しく説明されています。

誰かが私の問題に取り組みたいと思ったら、コードは以下にあり、グラフはdropbox.comのこのWordファイルにあります。はい、これは宿題の問題です。私の思考過程の何が悪いのかを理解したいので、これを投稿します。

これは私のRK4コードです：

これは私のオイラーコードです：

ルンゲクッタ法を使用してローレンツシステムを計算しようとしていますが、コードにエラーがある場所が見つかりません。それを実行すると、システムは静的なポイントに移動し、バタフライ（ローレンツアトラクター）を取得する必要があります。これは、ルンゲクッタ法のセクションの「for」ループにあるものだと思います。これが私のコードです

前もって感謝します

編集

Linuxで（より単純化された方法で）コードを再度記述しましたが、正常に動作します。Windowsのエディターには、コードをコンパイルしたときに多くの無限の値がスローされる何か（おそらくエンコーディング）があったようです。理由がわからないので、それに気付くのにかなりの時間がかかりました。ご協力いただきありがとうございます

これは私の最初の投稿であり、認めざるを得ません。プログラミングがひどいです。私はクラスのその男であり、彼の尻尾を切り落としていますが、他のクラスメートと同様にプログラミングを理解しているようには見えません。よろしくお願いします。以下で私の問題を説明しようと思います。

次のコード（コメントは削除されています）がありますが、実行すると、以下のような警告が表示されます。また、プログラムを実行すると、最初のユーザー入力値が許可されますが、突然、プログラムの最後にジャンプし、他の変数（変数「ベータ」など）の値を入力できなくなります。 ）。出力の画像（http://i.stack.imgur.com/yc3jq.jpg）があり、アルファを入力したことがわかりますが、プログラムは最後まで実行されます。何かご意見は？

どうもありがとうございました！-スペンサー

- - - - - - - - - - - - - - -コード - - - - - - - -

警告の例：

Does anyone know any c++/opengl sourcecode demos for 2D rigid body physics using runge kutta?

I want to build a physics engine but I need some reference code to understand better how others have implemented this.

C/c++ をやってからしばらく経ちましたが、次の ODE セットを解くために gsl ライブラリを使用して ODE ソルバーを書きたかったのです。

したがって、gsl 表記では my y[0]=u, y[1]==up となり、上記の RHS は f[0] と f[1] を定義します。これらの定義から、ヤコビアンと dfdr を計算できます (通常、それらの「時間」変数は「r」ではなく「t」と呼ばれます)。これを行う理由は、Mathematica の速度に問題があるためです。私は、ODE ソルバーに関するドキュメントの最後にある gsl サンプル コードを取り上げ、次のように私の問題に適用しようとしました。

これは数値を与えていますが、WorkingPrecision と PrecisionGoal が低い場合でも、Mathematica NDSolve が与える数値と同じではありません。私が行ったことに誤りはありますか？

この回答で私に与えられたアドバイスに従って、私は重力シミュレーターにルンゲクッタ積分器を実装しました。

ただし、ソーラーシステムの1年間をシミュレートした後でも、位置はcca 110 000 kmずれており、これは許容できません。

私の最初のデータはNASAのHORIZONSシステムによって提供されました。それを通して、私は特定の時点での惑星、冥王星、衛星、デイモス、フォボスの位置と速度のベクトルを取得しました。

これらのベクトルは3Dでしたが、惑星が太陽の周りのプレートに整列しているため、3次元を無視できると言われたので、そうしました。xy座標をファイルにコピーしただけです。

これは私の改善された更新方法のコードです：

アルゴリズムは次のとおりです。

1. システム内のすべてのボディの加速度を更新します
2. RK4（最初のステップ）
3. goto 1
4. RK4（秒）
5. goto 1
6. RK4（3番目）
7. goto 1
8. RK4（4番目）

RK4の実装で何かを台無しにしましたか？または、破損したデータから始めただけですか（重要なボディが少なすぎて、3次元を無視しています）？

これはどのように修正できますか？

私のデータなどの説明...

私の座標はすべて太陽を基準にしています（つまり、太陽は（0、0）にあります）。

110 000 kmシミュレーターで予測されたものからNASAによって与えられた地球のx座標を引くことによって、エラーが発生しました。

相対誤差はごくわずかに見えますが、それは単に、私のシミュレーションとNASAの両方で、地球が太陽から本当に遠く離れているためです。距離はまだ大きく、シミュレータが役に立たなくなります。

ベジエ曲線 があるとしましょう。パラメーターを一定の割合でB(u)インクリメントuすると、曲線に沿って一定速度の移動がu得られません。これは、パラメーターと曲線を評価して得られたポイントとの関係が線形ではないためです。

David Eberly の記事を読んで実装しました。パラメトリック曲線に沿って一定速度で移動する方法を説明します。

F(t)時間の値を入力として受け取る関数と、 timeでの速度の値を返すt速度関数があるとします。一定の割合で t パラメータを変化させて、曲線に沿って一定速度の移動を取得できます。sigmatB(F(t))

私が使用している記事のコアは、次の関数です。

提供された時間とシグマ関数uを使用して計算された曲線パラメーターを取得できます。tこれで、速度シグマが一定の場合に関数が正常に機能します。シグマが一様な加速を表している場合、それから間違った値が得られます。

これは、直線ベジエ曲線の例です。ここで、P0 と P1 はコントロール ポイントで、T0 T1 は接線です。曲線の定義:

time での曲線に沿った位置を知りたいとしましょうt = 3。速度が一定の場合:

および次のデータ:

位置を分析的に計算できます。

ベジエ スプラインと上記のアルゴリズムを使用して同じ方程式を解くと、次のようになりますn =5。

数値近似を考慮すると、値は非常に正確です (私はそれについて多くのテストを行いました。詳細は省略しますが、ベジエ曲線の実装は問題なく、曲線自体の長さはGaussian Quadratureを使用して非常に正確に計算されます)。

ここで、sigma を一様な加速関数として定義しようとすると、悪い結果が得られます。次のデータを検討してください。

線形運動方程式を使用して、粒子が P1 に到達する時間を計算できます。

私は加速度をt計算することができます：

シグマ関数を定義するためのすべてのデータがあります。

これを解析的に解決すると、時間後の粒子の速度は次のようになると予想されますt =3。

位置は次のようになります。

しかし、上記のアルゴリズムで計算すると、位置は次のようになります。

それは私が期待しているものとはかなり異なります。

長文失礼しました、どなたか読んでいただけると嬉しいです。

何か提案はありますか?何が欠けていますか?誰が私が間違っているのか教えてくれますか?

編集: 3D ベジエ曲線を試しています。次のように定義します。

および導関数: