MATLABコードヘルプ。後方オイラー法

StackOverflow https://stackoverflow.com/questions/2937183

質問

これが次のとおりです マトラブ/フリーメット コード私は解決しなければなりませんでした オード 数値的に使用して 後方オイラー 方法。しかし、結果は私の教科書の結果と矛盾しており、時には途方もなく矛盾しています。コードの何が問題になっていますか？

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
役に立ちましたか？

解決

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

フォワードオイラー： 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 間違った点で。

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

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

他のヒント

私が見つけることができる唯一の問題は、ライン：

n=(xfinal-xinit)/h

あるべきです：

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

負のステップ数を回避するため。ネガティブX方向に移動している場合は、関数に負のステップサイズを与えるようにしてください。

あなたの答えは、あなたがあなたの答えをどれほど粗く近似しているかのために、おそらく逸脱します。半acccurateの結果を得るには、デルタックスは非常に小さくなければならず、ステップサイズは非常に非常に小さくなければなりません。

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

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

見て 数値レシピ, 、具体的には第16章、通常の微分方程式の統合。オイラーの方法は問題があることが知られています。

eulerの方法が実際の使用には推奨されない理由はいくつかあります。

したがって、教科書がEulerの方法を使用していることを知らない限り、結果が一致するとは思わないでしょう。たとえそうであっても、同一の結果を得るには、同一のステップサイズを使用する必要があります。

あなたがいない限り 本当 あなたが自分で書いたオイラーの方法を介してオードを解決したいあなたはあなたが見るべきです 組み込みのODEソルバー.

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

Matlabが実装したソルバーは次のとおりです。

Ode45は、ドーマンドプランスペアである明示的なRunge-Kutta（4,5）式に基づいています。これはワンステップソルバーです。コンピューティングY（TN）では、直前の時点Y（TN-1）で解決策のみが必要です。一般に、ODE45は、ほとんどの問題の最初の試みとして適用するのに最適な機能です。

Ode23は、BogackiとShampineの明示的なRunge-Kutta（2,3）ペアの実装です。粗耐性や中程度の剛性の存在下では、ODE45よりも効率的かもしれません。 ODE45と同様に、ODE23はワンステップソルバーです。

ODE113は、Adams-Bashforth-Moulton Peceソルバーの可変順序です。厳しい許容範囲でODE45よりも効率的である可能性があり、ODEファイル関数が特に評価するのに費用がかかる場合。 ODE113はMultiStepソルバーです。通常、現在のソリューションを計算するためにいくつかの前の時点でソリューションが必要です。

上記のアルゴリズムは、非スティフシステムを解決することを目的としています。彼らが不当に遅いように見える場合は、以下の硬いソルバーのいずれかを使用してみてください。

ODE15Sは、数値分化式（NDFS）に基づく可変オーダーソルバーです。オプションで、通常は効率が低い後方分化式（Gearの方法としても知られているBDFS）を使用します。 ODE113と同様に、ODE15Sはマルチステップソルバーです。 ODE45が失敗した場合、または非常に非効率的である場合、ODE15Sを試してください。問題が硬いと思われます。

ODE23Sは、注文2の修正されたRosenBrock式に基づいています。1段階のソルバーであるため、粗耐性のODE15よりも効率的かもしれません。 ODE15が効果的でないいくつかの種類の硬い問題を解決できます。

ODE23Tは、「自由な」補間を使用した台形ルールの実装です。問題が適度に硬いだけで、数値減衰なしで解決策が必要な場合は、このソルバーを使用してください。 ODE23TはDAESを解くことができます。

ODE23TBは、TR-BDF2の実装であり、暗黙のRunge-Kuttaフォーミュラであり、第1段階の第1段階であり、2次の第2段階である第2段階です。構造により、同じ反復マトリックスが両方の段階の評価に使用されます。 ODE23Sと同様に、このソルバーは、粗耐性でODE15よりも効率的です。

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

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

コードは問題ありません。 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

ダミー関数をチェックしましたが、結果は有望でした。

ライセンス： CC-BY-SA と 帰属
所属していません StackOverflow
scroll top