2

ルンゲクッタ法（「RK4」）が常微分方程式の正確な解の0.01％以内になるために必要なステップ数を決定しようとしています。これをオイラー法と比較しています。どちらも、両対数プロット上で直線になるはずです。私のオイラーソリューションは正しいようですが、RKのカーブしたソリューションを取得しています。それらは同じコードに基づいているので、私は問題について完全に混乱しています。

編集：ウィキペディアのリンクを削除して申し訳ありません。私は新しいユーザーなので、複数のリンクを保持することはできません。ただし、どちらの方法も、私の実装と同様にWikipediaで詳しく説明されています。

誰かが私の問題に取り組みたいと思ったら、コードは以下にあり、グラフはdropbox.comのこのWordファイルにあります。はい、これは宿題の問題です。私の思考過程の何が悪いのかを理解したいので、これを投稿します。

f = @(x,y) x+y; %this is the eqn (the part after the @(t,y)

これは私のRK4コードです：

k1=@(x,y) h*f(x,y);
k2=@(x,y) h*f(x+1/2*h,y+1/2*k1(x,y));
k3=@(x,y) h*f(x+1/2*h,y+1/2*k2(x,y));
k4=@(x,y) h*f(x+h,y+k3(x,y));

clear y x exact i
x(1)=0;
y(1)=2;
xn=1;
x0=0;

exact=3*exp(xn)-xn-1; %exact solution at x=1
%# Evaluate RK4 with 1 step for x=0...1
N=1; %# number of steps
h=(xn-x0)/N; %# step size
i=1;
y(i+1)=y(i)+1/6*k1(x(i),y(i))+1/3*k2(x(i),y(i))+1/3*k3(x(i),y(i))+1/6*k4(x(i),y(i));
RK4(N)=y(i+1);  %# store result of RK4 in vector RK4(# of steps)
E_RK4(N)=-(RK4(N)-exact)/exact*100;%keep track of %error for each N
Nsteps_RK4(N)=N;


%# repeat for increasing N until error is less than 0.01%
while -(RK4(N)-exact)/exact > 0.0001
    i=1;
    N=N+1;
    h=(xn-x0)/N;
    for i=1:N
        y(i+1)=y(i)+1/6*k1(x(i),y(i))+1/3*k2(x(i),y(i))+1/3*k3(x(i),y(i))+1/6*k4(x(i),y(i));
        x(i+1)=x(i)+h;
    end
    RK4(N)=y(i+1);
    Nsteps_RK4(N)=N;
    E_RK4(N)=-(RK4(N)-exact)/exact*100; %# keep track of %error for each N
end

Nsteps_RK4(end);

これは私のオイラーコードです： 

%# Evaluate Euler with 1 step for x=0...1
clear y x i
x(1)=0;
y(1)=2;
N=1; %# number of steps
h=(xn-x0)/N; %# step size
i=1;
y(i+1)= y(i)+h*f(x(i),y(i));
Euler(N)=y(i+1); %# store result of Euler in vector Euler(# of steps)
E_Euler(N)=-(Euler(N)-exact)/exact*100;%# keep track of %error for each N
Nsteps_Euler(N)=N;
%# repeat for increasing N until error is less than 0.01%
while -(Euler(N)-exact)/exact > 0.0001
    i=1;
    N=N+1;
    h=(xn-x0)/N;
    for i=1:N
        y(i+1)= y(i)+h*f(x(i),y(i));   
        x(i+1)=x(i)+h;
    end
    Euler(N)=y(i+1);
    Nsteps_Euler(N)=N;
    E_Euler(N)=-(Euler(N)-exact)/exact*100; %# keep track of %error for each N
end
4

1 に答える 1

1

参照： http： //www.mathworks.com/help/techdoc/matlab_prog/f4-70115.html#f4-71621

式[...]の本体で指定された変数には、それらを使用する無名関数を作成するときに値が割り当てられている必要があります。構築時に、MATLABは、その関数の本体で指定された各変数の現在の値をキャプチャします。値がワークスペースで変更されたり、スコープ外になったりした場合でも、関数はこの値を変数に関連付け続けます。

ループ内で変更しても、 hin k1...の値は一定のままです。k4hwhile

h1つの解決策は、無名関数に追加することです。

k1=@(x,y,h) h*f(x,y);
于 2011-10-10T02:06:45.267 に答える