9

誰かが次の問題を手伝ってくれれば幸いです。次の ODE があります。

dr/dt = 4*exp(0.8*t) - 0.5*r   ,r(0)=2, t[0,1]       (1)

私は(1)を2つの異なる方法で解決しました。Runge-Kutta 法(4 次) とMatlabode45を使用。両方の結果を次の式で得られる分析解と比較しました。

r(t) = 4/1.3 (exp(0.8*t) - exp(-0.5*t)) + 2*exp(-0.5*t)

正確な解に対して各メソッドの絶対誤差をプロットすると、次のようになります。

RK メソッドの場合、私のコードは次のとおりです。

h=1/50;                                            
x = 0:h:1;                                        
y = zeros(1,length(x)); 
y(1) = 2;    
F_xy = @(t,r) 4.*exp(0.8*t) - 0.5*r;                   
for i=1:(length(x)-1)                              
    k_1 = F_xy(x(i),y(i));
    k_2 = F_xy(x(i)+0.5*h,y(i)+0.5*h*k_1);
    k_3 = F_xy((x(i)+0.5*h),(y(i)+0.5*h*k_2));
    k_4 = F_xy((x(i)+h),(y(i)+k_3*h));
    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  % main equation
end

ここに画像の説明を入力

そしてのためにode45

tspan = 0:1/50:1;
x0 = 2;
f = @(t,r) 4.*exp(0.8*t) - 0.5*r;
[tid, y_ode45] = ode45(f,tspan,x0);

ここに画像の説明を入力

私の質問は、なぜ使用すると振動するのode45ですか? (私は絶対誤差を指しています)。どちらの解も正確 ( 1e-9) ですがode45、この場合はどうなるでしょうか?

RK 法の絶対誤差を計算すると、見栄えが良くなるのはなぜですか?

4

2 に答える 2

22

RK4 関数は、実行中のステップよりもはるかに小さい固定ステップを実行していますode45。実際に表示されているのは、実際のステップ間のポイントを生成するために使用される多項式補間によるエラーですode45。これはしばしば「高密度出力」と呼ばれます ( Hairer & Ostermann 1990を参照)。

TSPAN3 つ以上の要素を持つベクトルを指定すると、 Matlab の ODE スイート ソルバーは固定ステップ サイズの出力を生成します。これは、彼らが実際に固定のステップ サイズを使用すること、または指定されたステップ サイズを使用することを意味するものではありませんTSPANode45実際に使用されたステップ サイズを確認し、構造体を出力して以下を使用することで、目的の固定ステップ サイズの出力を取得できますdeval

sol = ode45(f,tspan,x0);
diff(sol.x) % Actual step sizes used
y_ode45 = deval(sol,tspan);

の最初のステップの後、0.02ODE が単純である0.1ため、後続のステップで に収束することがわかります。これは、既定の最大ステップ サイズ制限 (積分間隔の 10 分の 1) と組み合わされた既定の許容誤差によって決まります。真のステップでエラーをプロットしましょう。

exactsol = @(t)(4/1.3)*(exp(0.8*t)-exp(-0.5*t))+2*exp(-0.5*t);
abs_err_ode45 = abs(exactsol(tspan)-y_ode45);
abs_err_ode45_true = abs(exactsol(sol.x)-sol.y);
abs_err_rk4 = abs(exactsol(tspan)-y);
figure;
plot(tspan,abs_err_ode45,'b',sol.x,abs_err_ode45_true,'k.',tspan,abs_err_rk4,'r--')
legend('ODE45','ODE45 (True Steps)','RK4',2)

エラープロット

ご覧のとおり、真のステップでのエラーは、RK4 のエラーよりもゆっくりと大きくなります (ode45実質的に RK4 よりも高次の方法であるため、これは予想されることです)。補間により、積分点の間で誤差が大きくなります。これを制限したい場合は、 で公差またはその他のオプションを調整する必要がありますodeset

ode45ステップの使用を強制したい場合は1/50、これを行うことができます (ODE は単純なので機能します)。

opts = odeset('MaxStep',1/50,'InitialStep',1/50);
sol = ode45(f,tspan,x0,opts);
diff(sol.x)
y_ode45 = deval(sol,tspan);

t = 10別の実験として、積分間隔を広げて多分に積分してみてください。エラーには多くの興味深い動作が見られます (ここでは相対エラーをプロットすると便利です)。これを説明できますか?ode45andを使用しodesetて、適切に動作する結果を生成できますか? 指数関数を広い間隔でアダプティブ ステップ法と統合することは困難でありode45、必ずしも最適なツールではありません。ただし、代替手段はありますが、プログラミングが必要になる場合があります。

于 2014-02-18T18:04:41.187 に答える
0

ode45 は rk4-rk5 に連結されています。個人的には ODE45 エラーの方がいいと思います。制限されたままであることに注意してください。エラーの大きさが大きくなりすぎると ode4 が修正され、1 サイクルあたりの最小エラーは約 1e-10 になります。rk4 は「逃げて」おり、それを止めるものは何もありません。

于 2014-02-18T16:51:36.210 に答える