問題タブ [odeint]

SciPy の integrate.odeint 関数を使用して、硬い ODE のシステムを統合しています。統合は簡単ではなく、時間がかかるため、対応する jacobian も使用しています。方程式を並べ替えることで、ヤコビアンを帯行列として定義できます。API ドキュメントに従って、muおよびmlパラメータで形状を定義したいと思います。残念ながら、ドキュメントは少しあいまいであるため、jacobian 関数を実装する方法を理解できませんでした。

odeintを呼び出す方法を確認するために、次の (ややばかげた) コードを使用しています。

完全な NxN ヤコビ行列を使用すると、統合は成功します。対角線のみを使用し、mu=0およびml=0を使用すると、積分も成功します。

バンド行列のユース ケースをテストするために、 mu=1とml= 1を使用して人工的な 3xN バンド ヤコビアンを作成しています。ここで、対角線からの導関数はすべてゼロです。これにより、ソルバーの奇妙な動作が発生します (非対角がゼロでない元の問題で見られるものと同様です)。

SciPyのintegrate.odeintでbanded jacobianオプションを使用する適切な方法は何ですか?

9x9 と 9x1 の 2 つの numpy 配列があります。離散時間点での微分方程式を解きたいのですが、ODEInt を機能させるのに問題があります。自分が正しいことをしているかどうかさえ確信が持てません。

Mathematica では、方程式は次のようになります。

Ab は 9x9 マトリックス、A0 は 9x1 マトリックス (初期) です。ここで、私は時間を解決し、人生は良いです.

Python の実装では、間違った答えを与える次のコードがあります。

出力は次のとおりです。

しかし、正しい答えは次のとおりです。

誰かが私が間違っているかもしれないことを指摘できますか?

3 つの複素行列と、これらの行列を結合した連立微分方程式があるとします。

基本的に2 つの質問があります。

1. コードを減らす方法は? reduce とは、すべてのコンポーネントを個別に書き留めるのではなく、ODE のシステムを解きながら行列を処理することで、これを行う方法はありますか?

2. t (私のコードの最後の 2 行)に関して行列の要素をプロットする代わりに、どうすれば固有値 (絶対値) (行列 A の固有値の abs を t の関数として) をプロットできますか?

私は odeint と推力の上に分析ツールを構築しています。このツールは、多数の初期条件の問題を解決します。私は odeint と throws のチュートリアル / デモに従って成功しました。これまで私が扱ってきたシステムには 5 つ未満の状態変数しかありませんでしたが、より多くの状態変数を処理できるようにシステムを拡張する必要があります。

チュートリアルで述べたように、thrust のタプルは最大 10 個のアイテムしか持てないという小さな問題があります。この問題は、ネストされた zip_iterators を使用することで解決できます。

チュートリアルの1つの下部近くから引用：

「Thrust のタプルの最大アリティが 10 であるという小さな問題があります。しかし、zip イテレータを詰め込んだ zip イテレータを作成できるので、これは小さな問題にすぎません。そのため、最上位の zip イテレータには、状態用の 1 つの zip イテレータと 1 つの通常のイテレータが含まれています。パラメータには 1 つの zip イテレータ、派生物には 1 つの zip イテレータを使用します。"

このソリューションを実装しましたが、機能しますが、残念ながら速度が大幅に低下します。同じシステム (8192 個の初期条件を同時に解く 2 変数システム) の場合、元のシンプルだが 5 つ以上の変数に拡張できないソリューションは 1 秒強で実行されます。

一方、より複雑ではあるが拡張可能なネストされたソリューションには、最大 2000 倍の時間がかかります!

2 つのプログラムの唯一の違いは、

1. functor の内部でoperator()、状態変数と導関数を参照し、

2. ファンクターのコンストラクター、特にfor_eachコマンド内で、zip_iterators とタプルを作成します。

以下に、各プログラムについて、これらのセクションからの抜粋を示します。この種の速度低下は壊滅的なので、ここで何か間違ったことをしているといいのですが! どんな助けでも大歓迎です。

「SIMPLE」コードの抜粋 (ネストされていないイテレータ)

「EXTENSIBLE」コードの抜粋 (ネストされたイテレータ)

odeint C++ ライブラリrunge_kutta4のメソッドを使用したいと考えています。Matlabで問題を解決しました。を初期値,で解決するための Matlab の私の次のコードは次のとおりです。x'' = -x - g*x'x1 = 1x2 = 0

main.m

ODESolver.m

odeint ライブラリをインストールしました。使用するための私のコードrunge_kutta4は次のとおりです

その結果が次の写真です

私の質問は、なぜ結果が異なるのですか? 私の C++ コードに何か問題がありますか?

これらは両方の方法の最初の値です

アップデート：

問題は、C++ コードに初期値を含めるのを忘れたことです。それに気づいてくれた@horchlerに感謝します。適切な値を含め、runge_kutta_dopri5@horchler が提案したように使用すると、結果は次のようになります。

これらの変更を反映するようにコードを更新しました。

私はPythonが初めてで、scipyの組み込みodeintの代わりにPythonでmatlab ode15sを使用して模倣したいと考えています。

コードはもともと次のように書かれています。

ここで、PSP は次のように定義されます。

scipy.integrate.odeint() の説明を読むと、この関数は次の形式で引数を取ります。

これは、func=PSP, y0=Rphi, t=t, args=(b,k,f) So Rphi が PSP 関数に入り、積分されて yy になって出てくることを意味し、この関数は t のすべての要素に対してそれを繰り返します。

今、それを matlab の ode15s を模倣するものに変換したいと考えています。他のいくつかのトレッドを読んで、私はそれを使用できることがわかりました

問題は、元の引数をこのインテグレータに渡すにはどうすればよいかということです。

おそらく次のようになると思います。

f が私の PSP 関数であることはわかっていますが、y0 は同じです: Rphi, ここに私の質問があります:

1. t0 の初期値は何ですか? t[0] だけですか?
2. 変数 (b,k,f) を関数 f=PSP に渡すにはどうすればよいですか?
3. これを ode15s と呼ぶ場合、t のベクトル サイズを積分して yy の最終値を収集するにはどうすればよいですか?

どんな助けでも大歓迎です。ありがとうございました。