私は次のpython MWEを持っています（コードは以下で説明されています）

ご覧のとおりbase_equations、本質的にガウス パルスである一連の方程式 ( ) を統合しています。odeintパルスの 2 つの中心点 (30 と 60) についてこれらの方程式を数値的に解くために使用します。

最初の中心点 (t=30) では、式 y1 により期待される動作が得られます。つまり、パルスが表示されます。

2 番目の中心点 (t=60) では、式 y2 によって予期しない動作が発生します。パルスがまったく表示されません!

働くことと働かないことの切り替えは、47 歳から 48 歳の間に起こります。

グラフィカルな出力は次のとおりです。予想される出力は、ライン y2a と y2b が 60 付近で劇的な変化を示すことですが、そうではありません。

何が起こっているのでしょうか？

私はHaskellを学んでいて、Cでできる限り速くコードを書こうとしています。この演習では、単純な1次元物理システム用のオイラー積分器を作成しています。

• CコードはGCC4.5.4およびでコンパイルされてい-O3ます。1.166秒で実行されます。
• HaskellコードはGHC7.4.1とでコンパイルされてい-O3ます。21.3秒で実行されます。
• Haskellをでコンパイルすると、 4.022秒-O3 -fllvmで実行されます。

それで、Haskellコードを最適化するための何かが欠けていますか？

PS：私は次の引数を使用しました：1e-8 5。

Cコード：

Haskellコード：

私は C++ で結合された ODE のセットを解こうとしていますが、odeint (ここで入手可能) について多くの良いことを読んだので試してみたいと思いました。ただし、パッケージに含まれている例をコンパイルしています。状況は次のとおりです。

これまでのところ、パッケージをダウンロードして解凍しただけです。Windows XP で Visual C++ 6 を実行しており、サンプル "simple1d" をコンパイルしています。「インクルード ファイルを開けません: 'boost/config.hpp': そのようなファイルまたはディレクトリはありません」というエラー メッセージが表示されます。

私が欠けているものは私には明らかではありません。これを機能させるには完全なブースト ライブラリが必要ですか? ODEINT は「スタンドアロン」パッケージであるという印象を受けました。

ode45などのODEソルバーを使用しても、呼び出された関数内のパラメーターの値を「変更」することはできますか？

たとえば、次の関数を使用する場合：

このスクリプトで：

すべてが期待どおりに機能します。ただし、パラメーター値を簡単に変更する方法が必要ですが、それでも1つの関数と1つのスクリプトだけでODEソルバーの複数の反復を実行します。ただし、ODEソルバーに用語を追加するだけでは不十分なようです。たとえば、次のようになります。

このスクリプトで：

動作しません。私はこれがうまくいかないことを知っていると思いますが、私は解決策を思い付くことができませんでした。また、関数内でifステートメントを検討し、使用するパラメーターのセットをいくつかハードコーディングしました（たとえば、P == 1の場合はセット1を使用し、P == 2の場合はセット2を使用します）が、これも次のようには機能しませんでした。 ODEで使用するセットをどこで呼び出すかはわかりません。パラメータ値を変更できるようにしながら、ODEソルバーで1つの関数と1つのスクリプトを使用する方法に関するヒントや提案をいただければ幸いです。

ありがとう、マイク

の (簡単な) ドキュメントにscipy.integrate.odeは、2 つのメソッド (dopri5およびdop853) にはステップサイズ制御と高密度出力があると記載されています。例とコード自体を見ると、インテグレーターから出力を取得するための非常に単純な方法しかわかりません。つまり、固定された dt だけ積分器を進め、その時点で関数値を取得し、繰り返すように見えます。

私の問題はかなり可変的なタイムスケールを持っているので、必要な許容範囲を達成するために評価する必要がある時間ステップで値を取得したいと思います。つまり、早い段階で物事がゆっくりと変化しているため、出力の時間ステップが大きくなる可能性があります。しかし、物事が面白くなるにつれて、出力のタイム ステップを小さくする必要があります。実際には、等間隔で高密度の出力は必要ありません。適応関数が使用する時間ステップが必要なだけです。

編集：密な出力

関連する概念 (ほぼ反対) は「高密度出力」です。これにより、実行されるステップは、ステッパーが実行するのと同じくらい大きくなりますが、関数の値は (通常、ステッパーの精度に匹敵する精度で) 補間されます。あなたがしたい。基礎となる fortranscipy.integrate.odeは明らかにこれが可能ですodeが、インターフェースがありません。 odeint一方、は別のコードに基づいており、明らかに高密度の出力を行います。(右側が呼び出されるたびに出力して、それがいつ発生するかを確認し、出力時間とは関係がないことを確認できます。)

したがって、事前に必要な出力時間ステップを決定できる限り、適応性を利用することができます。残念ながら、私のお気に入りのシステムでは、統合を実行するまで、時間の関数としてのおおよそのタイムスケールが何であるかさえわかりません。したがって、積分器を 1 ステップ実行するというアイデアと、この高密度出力の概念を組み合わせる必要があります。

EDIT 2:再び密な出力

どうやら、scipy 1.0.0 では、新しいインターフェイスを介した高密度出力のサポートが導入されました。特に、キーワードとしてのから離れscipy.integrate.odeintたり近づいたりすることを推奨しています。に設定すると、返されるオブジェクトには、時間の配列で呼び出すことができる属性が含まれます。この属性は、それらの時間に統合された関数の値を返します。それでもこの質問の問題は解決しませんが、多くの場合に役立ちます。scipy.integrate.solve_ivpdense_outputTruesol

おそらくGSLライブラリを使用してCで書かれたODEソルバーは、Mathematica 8.0 NDSolveと比較して大幅な速度の利点がありますか？精度の点でどのように公平でしょうか？

私の理解では、コンパイルされたコードは原則としてより高速である可能性がありますが、最近のNDSolveは、すでに何らかの形でコンパイルされたコード自体を大量に使用していますか？

また、MathLinkやMathematicaのコンパイル関数などを使用して、ODEの解決を高速化するためのオプションはありますか？

私はいくつかの頌歌を含む次の関数を持っています:

以下にいくつかの結果を示します。

州 Y は異なります。理由を説明してもらえますか?

しきい値を正しく設定していないと思います。その方法はありますか？

ありがとう！

したがって、初期条件 x(0)= 0 および v(0) = x'(0) = v_o = 1 で x''(t) = -x(t)^p を解く必要があります。パラメーター p の値は 1 です。

これは私が持っているものです:

これは私が受け取るエラーメッセージです：

ode_oscillation(1) odeearguments の使用エラー (91 行目) ODE_OSCILLATION/OSCILATION_EQUATION は列ベクトルを返さなければなりません。

ode45 のエラー (114 行目) [neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

ode_oscillation のエラー (17 行目) [t,velocity,x] = ode45(@oscilation_equation, tspan,1);

単純なifループを使用して、odeスクリプト内のパラメーター値を変更しています。これは私が書いた同じ問題を示すスクリプトの例です。したがって、最初に機能するバージョンは次のとおりです。

これは、を使用して実行できます

と視覚化

わかりました、それは良いバージョンです。今、あなたがするのがtをに変更することだけなら

全部が爆破します。グラフを平均化するのはMATLABだと思うかもしれませんが、それは

コードが変更されていないため、どちらの場合も答えは同じであるはずです。しかし、この2番目のバージョンは0を出力しますが、これは間違っています。

一体何が起こっているのですか？誰かアイデアがありますか？

助けてくれてありがとう

次のようなordinary differential equationインタイムドメインがあります。

どこI = A*t/tau*exp^(1-t/tau)

周波数ドメイン:

j複素数であることsqrt(-1)

したがって、高速フーリエ変換( )u(t)を使用して周波数領域に入り、 を使用し て戻ることができます。fftifft

コード:

しかし、これu(t)を微分方程式の数値積分やその解析形式などの他の方法と比較すると、正しくありません。私は自分の間違いがどこにあるかを理解しようとしましたが、うまくいきませんでした。

啓発してください。