0

これは、後方オイラー法を使用してODEを数値的に解くために取得したMATLAB / FreeMatコードです。しかし、結果は私の教科書の結果と一致せず、時にはばかげて一貫性がないことさえあります。コードの何が問題になっていますか？

function [x,y] = backEuler(f,xinit,yinit,xfinal,h)

    %f - this is your y prime
    %xinit - initial X
    %yinit - initial Y
    %xfinal - final X
    %h - step size

    n = (xfinal-xinit)/h; %Calculate steps

    %Inititialize arrays...
    %The first elements take xinit and yinit corespondigly, the rest fill with 0s.
    x = [xinit zeros(1,n)];
    y = [yinit zeros(1,n)];

    %Numeric routine
    for i = 1:n
        x(i+1) = x(i)+h;
        ynew = y(i)+h*(f(x(i),y(i)));
        y(i+1) = y(i)+h*f(x(i+1),ynew);
    end
end
4

6 に答える 6

5

あなたのメソッドは新しい種類のメソッドです。後方オイラーでも前方オイラーでもありません。:-)

フォワードオイラー：y1 = y0 + h*f(x0,y0)

後方オイラーsolve in y1: y1 - h*f(x1,y1) = y0

あなたの方法：y1 = y0 +h*f(x0,x0+h*f(x0,y0))

あなたの方法は後方オイラーではありません。

で解くのではなく、フォワードオイラー法でy1推定するだけです。私はあなたの方法の分析を追求したくありませんが、あなたが間違った点でy1関数を評価するので、フォワードオイラーと比較しても、それは確かにうまく機能しないと思います。f

これが私が考えることができるあなたの方法に最も近い方法であり、明示的でもあり、はるかに良い結果をもたらすはずです。それはホイン法です：

y1 = y0 + h/2*(f(x0,y0) + f(x1,x0+h*f(x0,y0)))

于 2010-05-30T07:29:56.207 に答える
2

私が見つけることができる唯一の問題は、次の行です。

n=(xfinal-xinit)/h

する必要があります：

n = abs((xfinal-xinit)/h)

負の歩数を避けるため。負のx方向に移動する場合は、関数に負のステップサイズを指定してください。

あなたの答えは、あなたがあなたの答えをどれほど粗く近似しているのかという理由でおそらく逸脱しています。半正確な結果を得るには、deltaXを非常に小さくする必要があり、ステップサイズを非常に小さくする必要があります。

PS。これは「逆オイラー法」ではなく、通常の古いオイラー法です。

これが宿題の場合は、タグを付けてください。

于 2010-05-30T01:20:34.343 に答える
2

数値レシピ、特に第16章、常微分方程式の積分をご覧ください。オイラー法には問題があることが知られています。

オイラー法が実際の使用に推奨されない理由はいくつかあります。その中には、（i）同等のステップサイズで実行される他のより洗練された方法と比較した場合、この方法はあまり正確ではなく、（ii）どちらも非常に安定していません。

したがって、教科書がオイラー法を使用していることを知らない限り、結果が一致するとは思われません。たとえそうだとしても、同じ結果を得るには、おそらく同じステップサイズを使用する必要があります。

于 2010-05-30T01:22:36.000 に答える
1

自分で作成したオイラー法を使用してODEを本当に解きたい場合を除いて、組み込みのODEソルバーを確認する必要があります。

x(i)補足：次のようにループ内に作成する必要はありませんx(i+1) = x(i)+h;。代わりに、単に書くことができますx = xinit:h:xfinal;。n = round(xfinal-xinit)/h);また、警告を避けるために書くこともできます。

これがMATLABによって実装されたソルバーです。

ode45は、明示的なルンゲクッタ（4,5）式、ドルマンプリンスペアに基づいています。これはワンステップソルバーです。y（tn）の計算では、直前の時点y（tn-1）の解のみが必要です。一般に、ode45は、ほとんどの問題の最初の試行として適用するのに最適な関数です。

ode23は、BogackiとShampineの明示的なルンゲクッタ（2,3）ペアの実装です。粗い公差で、適度な剛性が存在する場合は、ode45よりも効率的である可能性があります。ode45と同様に、ode23はワンステップソルバーです。

ode113は、可変次数のAdams-Bashforth-MoultonPECEソルバーです。厳しい公差で、ODEファイル関数の評価に特に費用がかかる場合は、ode45よりも効率的です。ode113はマルチステップソルバーです。通常、現在の解を計算するには、先行するいくつかの時点での解が必要です。

上記のアルゴリズムは、非剛性システムを解決することを目的としています。過度に遅いと思われる場合は、以下のスティッフソルバーのいずれかを使用してみてください。

ode15sは、数値微分公式（NDF）に基づく可変次数ソルバーです。オプションで、通常は効率が低い後退微分法（BDF、Gearの方法とも呼ばれます）を使用します。ode113と同様に、ode15sはマルチステップソルバーです。ode45が失敗するか、非常に非効率的で、問題が難しいと思われる場合、または微分代数問題を解く場合は、ode15sを試してください。

ode23sは、次数2の修正されたRosenbrock式に基づいています。これはワンステップソルバーであるため、大まかな許容誤差でode15sよりも効率的である可能性があります。ode15sが効果的でないある種の難しい問題を解決することができます。

ode23tは、「フリー」補間を使用した台形公式の実装です。問題が適度に硬く、数値減衰のない解が必要な場合は、このソルバーを使用してください。ode23tはDAEを解決できます。

ode23tbは、TR-BDF2の実装です。これは、台形公式ステップである第1ステージと、2次の後退微分式である第2ステージを持つ暗黙のルンゲクッタ式です。構造上、両方の段階の評価に同じ反復行列が使用されます。ode23sと同様に、このソルバーは、粗い許容誤差でode15sよりも効率的である可能性があります。

于 2010-05-30T01:40:45.417 に答える
0

このコードは機能すると思います。これを試して。

for i =1:n
    t(i +1)=t(i )+dt;
    y(i+1)=solve('y(i+1)=y(i)+dt*f(t(i+1),y(i+1)');
    end
于 2016-03-15T04:37:30.827 に答える
0

コードは問題ありません。forループ内に別のループを追加する必要があります。一貫性のレベルを確認します。

if abs((y(i+1) - ynew)/ynew) > 0.0000000001 
    ynew = y(i+1);
    y(i+1) = y(i)+h*f(x(i+1),ynew);
end

ダミー関数をチェックしたところ、期待できる結果が得られました。

于 2016-05-17T15:38:50.670 に答える