ルンゲクッタ法の入力引数の不足について。

function xy = rk4(fun,x0,xe,h,y0)%関数、初期値、終端値、分割数、y初期値
x  = linspace(x0,xe,h);%%h個の点、間隔は(xe-x0)/h
f = @fun(x);
y = zeros(1,h);
dx = (xe-x0)/h;
y(1) = y0;

Y=feval(f,x);
    
for k=2:h
    k1 = h*Y(k);
    k2 = h*feval(f,x(k)+k1./2);
    k3 = h*feval(f,x(k)+k2./2);         %%(Y(k)+k2.*dx./2)./(x(k)+dx./2);
    k4 = h*feval(f,x(k)+k3);
    y(k)=y(k-1)+(k1+2*k2+2*k3+k4).*dx/6;
end

cosを関数として入力した場合、入力引数が不足しています。

ある関数を入力すれば、その関数によってルンゲクッタ法を再現できるようにしたいのですが、うまく作ることができません。よろしくお願いします。

Answer 1

[michio]

関数自体を入力引数とする場合には、関数ハンドルとして入力するのがオススメです。

https://jp.mathworks.com/help/matlab/matlab_prog/creating-a-function-handle.html

例えば

xy = rk4(@cos, x0, xe, h, y0)

という風に対象となる関数の頭に @ を付けて入力します。

その場合、rk4 の関数内では fun(x) と fun という名前の関数として使用できます。（以下の「変更したところ」を確認ください）

function xy = rk4(fun,x0,xe,h,y0)%関数、初期値、終端値、分割数、y初期値
x  = linspace(x0,xe,h);%%h個の点、間隔は(xe-x0)/h
f = fun(x); % 変更したところ
y = zeros(1,h);
dx = (xe-x0)/h;
y(1) = y0;

Y=feval(f,x);
    
for k=2:h
    k1 = h*Y(k);
    k2 = h*feval(f,x(k)+k1./2);
    k3 = h*feval(f,x(k)+k2./2);         %%(Y(k)+k2.*dx./2)./(x(k)+dx./2);
    k4 = h*feval(f,x(k)+k3);
    y(k)=y(k-1)+(k1+2*k2+2*k3+k4).*dx/6;
end

[yamamoto yosuke]

返信ありがとうございます。 この変更をさせていただいて、入力をそれぞれ @cos,0,2＊pi,100,1としたところ、 fevalのエラーで、評価対象の関数は、string スカラー、文字ベクトルまたはfunction_handle オブジェクトであらわさなければなりません。となります。

[michio]

確かにその先の feval まで見ていませんでした。

せっかくですので、私がイメージしていたルンゲクッタを下記に記載しますので、参考になりましたら幸いです。 fun という入力引数は関数ハンドルを使用する変更により、それ自体を既に「関数」として使えるものなので、fun を特定の x で評価するにあたって feval は不要です。

y = rk4_(@cos,0,2*pi,100,1);
plot(y)

と実行してみてください。

function y = rk4_(fun,x0,xe,h,y0)%関数、初期値、終端値、分割数、y初期値
x  = linspace(x0,xe,h);%%h個の点、間隔は(xe-x0)/h
y = zeros(1,h);
dx = (xe-x0)/h;
y(1) = y0;
for k=2:h
    k1 = fun(x(k));
    k2 = fun(x(k)+dx*k1./2);
    k3 = fun(x(k)+dx*k2./2);         %%(Y(k)+k2.*dx./2)./(x(k)+dx./2);
    k4 = fun(x(k)+dx*k3);
    y(k)=y(k-1)+(k1+2*k2+2*k3+k4).*dx/6;
end

[yamamoto yosuke]

この変更でエラーなしに実行することを確認しました。結局、fevalを使わなければよかったのですね。参考になりました。ありがとうございました。

[michio]

よかったです！

もともと提示いただいたルンゲクッタの各ステップで、h がかけられていた部分が気になりました。私が記載したコードでは削除していますので、正しい計算になっているかどうかご確認くださいませ。よろしくお願いします。

[yamamoto yosuke]

そうですね、hについては、これから考えてみます。