問題タブ [odeint]

ODE in pythonを解決して結果をグラフ化するスクリプトを作成しようとしています。このタスクには scipy.integrate.odeint を使用しています。簡単なチュートリアルに従い、コードを修正して、解決したい ODE で動作するようにしました

このコードは、-x のような単純な関数に対しては完全に実行されますが、問題は、私が使用している関数に対しては機能しないことです。このメソッドからプロットを取得できますが、使用しても同じプロットが得られます

y=(1/(p*V*C)) * (A*e*b*np.power(T, 4) - (A*e*b*np.power(x, 4)) - ((1.32*A)/(D**0.25))*(np.power((x-T), 1.25)))

また

y=(1/(p*V*C)) * (A*e*b*pow(T, 4) - (A*e*b*pow(x, 4)) - ((1.32*A)/(D**0.25))*(pow((x-T), 1.25)))

また

y=(1/(p*V*C)) * (A*e*b*np.power(T, 4) - (A*e*b*np.power(x, 4))

これでも同じプロットが得られると思います

y=(1/(p*V*C)) * ((-A*e*b*np.power(x, 4))

また、e の値は 0 から 1 の間である必要がありますが、指数関数的減衰のように見えるものを得るには巨大な数を使用する必要があります。私は基本的にこれをやろうとしています: sfu.ca/~rca10/rants/convection.pdf. 問題の ODE は 4 ページの上部にあります。リンク内の人が通常の放射率を使用できる場合、すべてに対して同じ ODE を同じ値で解いていたとしても、私は使用できません (私はまったく同じラボを行っています)。

私は一体何を間違っているのですか？

そのため、3 つの ODE のシステムを解こうとしており、ODEint を使用してそれらを解く次のコードを開発しました。しかし、実行すると、ODEint で連立方程式の関数を呼び出すときに問題が発生します。

そして、私の出力では次のようになります。

y0 と t の順序を変更しても、同じエラーが発生します。私はエラーを理解していませんし、問題へのガイドもありません。どんな助けでも大歓迎です。

最初のモジュールは次の 3 つです。

私の2番目のもの：

そして私の3番目のもの：

3番目のメインモジュールからプログラムを実行しようとすると、次のエラーが発生します: in func a_x = -(220**2)*x/N - (G * MBH * x)/NT TypeError: can'tシーケンスに 'tuple' 型の非 int を掛ける

だから私は、プログラムが a_x がタプルであると考えていると思います。

このコードは、以前は a_x の (G*MBH*x)/NT 部分がなくても完全に機能していました。したがって、プログラムが a_x をシーケンスと見なす理由を理解する必要があると思います。

したがって、基本的には odeint を使用して自分が持っている頌歌を解き、それを別の方程式の変数として使用したいと考えています (最終的には実際に私の解である方程式)。

したがって、ここで最初の頌歌を解きます (私のコードは実際には ODE の行列を解いていますが、同じ問題が引き続き適用されます)。

私の f() 関数は次のように定義されています。

それで、すべての良いものの後、私は解決策を示すためにプロットできる配列または数値、ポイントを取得します。しかし、私が持っている質問は、この解決策を使用する必要があるということです（簡単にするために、ソリューションのマトリックスではありません）私が実際に見つけようとしている機能にプラグインします：

ここで、bernoulli() はその次数のベルヌーイ数を返す単なる関数であり、adjointOp() は指定された次数の随伴演算子 (交換関係) です。

OmegのソリューションがInterpolationFunctionであったMathematicaでこれを最初に行いましたが、それは単に変数として機能し、乗算などを行いました。ただし、Pythonでこの状況を処理する方法がわかりません。

アップデート：

この例を他のシステムで実行しました。Intel i7-3630QM、Intel HD4000、および Radeon HD 7630M では、すべての結果は同じです。i7-4700MQ / 4800MQ では、OpenCL または 64 ビット gcc を 32 ビット gcc から使用すると、CPU の結果が異なります。これは、デフォルトで SSE を使用する 64 ビット gcc および OpenCl と、387 math を使用する 32 ビット gcc の結果です。-mfpmath=387 が設定されている場合、少なくとも 64 ビット gcc は同じ結果を生成します。そのため、さらに多くの情報を読み、x86 浮動小数点を試してみる必要があります。回答ありがとうございます。

「CUDA と OpenCL のプログラミング: 最新の C++ ライブラリを使用したケーススタディ」の Lorenz システムの例を、それぞれ異なる OpenCL デバイス上の 10 個のシステムに対して実行しましたが、異なる結果が得られています。

1. Quadro K1100M (NVIDIA CUDA)

   r => xyz
   0.100000 => -0.000000 -0.000000 0.000000
   5.644444 => -3.519254 -3.5192504.644452 11.188890
   => 5.212534 5.212530 10.188904 16.73334
   = 6.4734 = 6.4734 = 6.4734 = 6.4734 = 6.473433737373737373737373737373734367373737347334>
   

   22.277779 => 3.178553 2.579687 17.946903
   27.822224 => 5.008720 7.753564 16.377680
   33.366669 => -13.381100 -15.252210 36.107887
   38.911114 => 4.256534 6.813675 ​​23.838787
   44.455555 => -11.083726 0.691549 53.632290
   50.000000 => -8.624105 -15.728293 32.516193
   

2. インテル(R) HD グラフィックス 4600 (インテル(R) OpenCL)

   r => xyz
   0.100000 => -0.000000 -0.000000 0.000000
   5.644444 => -3.519253 -3.519250 4.644451
   11.1888890 => 5.212531 5.212538 10.188888888888888888888890
   16.733334 = 15.47326> 15.47326> 15.47326
   

   22.277779 => 7.246771 7.398651 20.735369
   27.822224 => -6.295782 -10.615027 14.646572
   33.366669 => -4.132523 -7.773201 14.292910
   38.911114 => 14.183139 19.582197 37.943520
   44.455555 => -3.129006 7.564254 45.736408
   50.000000 => -9.146419 -17.006729 32.976696
   

3. Intel(R) Core(TM) i7-4800MQ CPU @ 2.70GHz (Intel(R) OpenCL)

   r => xyz
   0.100000 => -0.000000 -0.000000 0.000000
   5.644444 => -3.519254 -3.519251 4.644453
   11.188890 => 5.212513 5.212507 10.188900
   16.733334 = 6.4734 = 6.4734 =>
   

   22.277779 => -8.295195 -8.198518 22.271002
   27.822224 => -4.329878 -4.022876 22.573458
   33.366669 => 9.702943 3.997370 38.659538
   38.911114 => 16.105495 14.401397 48.537579
   44.455555 => -12.551083 -9.239071 49.378693
   50.000000 => 7.377638 3.447747 47.542763
   

ご覧のとおり、3 つのデバイスは R=16.733334 までの値で一致し、その後発散し始めます。

VexCL なしで odeint を使用して同じ領域を実行し、CPU 実行で OpenCL の結果に近い結果を得ました。

バニラ風味:

サンプル コードは次の場所にあります: https://github.com/ddemidov/gpgpu_with_modern_cpp/blob/master/src/lorenz_ensemble/vexcl_lorenz_ensemble.cpp

ここで何を見ているのかわからないのですか？CPU の結果は互いに非常に近いため、GPU の問題のように見えますが、私は OpenCL の初心者なので、この根本的な原因を見つける方法についていくつかの指針が必要です。

次の連立方程式の解を実装しました

最初は Mathematica で、その後 Python で。

Mathematica での私のコード:

そこから、次のプロットを取得します: Plot1 (403 Forbidden メッセージが表示された場合は、URL フィールド内で Enter キーを押してください)

後で、同じものを python にコーディングしました。

そして、これは私が得るプロットです: Plot2 (403 Forbidden メッセージが表示された場合は、URL フィールド内で Enter キーを押してください)

Mathematica の解をより小さなフィールドで再度プロットすると、次のようになります。

Python ソリューションと同様の結果が得られます。軸だけがずれます。

プロットに大きな違いがあるのはなぜですか？私は何を間違っていますか？

特に f1 が計算される場所で、モデルの Python 実装が間違っていると思われます。あるいは、この場合のように、パラメトリック方程式をプロットするのに plot() 関数がまったく役に立たないかもしれません。

ありがとう。

追伸：テキスト内の画像を平手打ちしないことであなたの人生を難しくして申し訳ありません。まだ評判が足りない。

私は、化学反応ネットワークを統合するために c++ ブースト/デイントを使用しています。一部の反応を排除する場合を除いて、すべて正常に機能します。次に、出力に -nan が含まれることがあります。これはその出力のごく一部です。最初の列は時間で、他の列は種の濃度です。

2 番目の種はそれ自体にのみ依存し、7 番目の種は定数です。

だから私の質問は: odeint はいつ出力として -nan を与えるのですか? どこかに 0/0 除算があるかもしれませんが、すべての方程式には積しか含まれていません。

dxdt[5] = - c[2] * x[5] + c[3] * x[9] * x[3];

アドバイスをありがとうございました。