0

Runge-Kutta4 を使用して、ODE ソリューション用に次のコードを作成しました。

function y=myODE(f,y0,x0,h,x_final)
x=x0;
y=y0;
while x<x_final
  h=min(h,x_final-x);
  k1=f(x,y);
  k2=f(x+h/2,y+h*k1/2);
  k3=f(x+h/2,y+h*k2/2);
  k4=f(x+h,y+h*k3);
  y=y+(h/6)*(k1+2*k2+2*k3+k4)
  x=x+h;
end

fは関数y' = f(x,y)、y0は初期値、x0は関数の開始位置、hサブインターバル、x_finalは関数の停止位置です。

コードを試してみたところ、ODE が正しく解決されましたが、サブインターバルx0を使用して、間隔の xy 軸上に関数をプロットしたいと考えています。を使用してプロットしようとすると、空のグラフしか表示されません。プロットするために複数にバインドする必要があることを（推測して）理解していますが、コードをあまり変更せずにそれを行うにはどうすればよいですか?x_finalhplot(x0:h:x_final,y)yx

ygiven y0、 interval x0to x_final、 givenのグラフをプロットするにはどうすればよいhですか?

MATLAB は初めてなので、できる限りのサポートをお願いします!

編集:私のコードの目的を明確にするため。

ソリューションとグラフ作成の両方に、この ODE ソルバーが必要です。yonの値を とh比較して切り捨て誤差を研究し、 のグラフを異なるで比較し2*hて Runge-Kutta4 の安定性を研究することになっています。 yh

4

3 に答える 3

1

これはコードのあまりスマートなリファクタリングではありません (ODE のステップ数によっては、解決が遅くなり、グラフィックスが殺されるため) ですが、眠いのでハックします:

    function y=myODE(f,y0,x0,h,x_final)
            hold(axes('Parent',figure),'on');
            x=x0;
            y=y0;
            plot(x,y, 'o');
            while x<x_final
                    h=min(h,x_final-x);
                    k1=f(x,y);
                    k2=f(x+h/2,y+h*k1/2);
                    k3=f(x+h/2,y+h*k2/2);
                    k4=f(x+h,y+h*k3);
                    y=y+(h/6)*(k1+2*k2+2*k3+k4);
                    x=x+h;
                    plot(x,y,'o');
            end;
    end

たぶん明日、この答えを適切なものに書き直すかもしれませんが、今のところ、おやすみなさい！:-)

于 2016-10-06T01:44:16.703 に答える
0 
function y=myode(f,y0,x0,h,x_final)
x=x0;
y=y0;
plot(x0,y0,'.')
hold on
while x<x_final
  h=min(h,x_final-x);
  k1=f(x,y);
  k2=f(x+h/2,y+h*k1/2);
  k3=f(x+h/2,y+h*k2/2);
  k4=f(x+h,y+h*k3);
  y=y+(h/6)*(k1+2*k2+2*k3+k4);
  x=x+h;
  plot(x,y,'.')
  disp([x,y])
end

コメントボックスでは、固定コードを「コード形式」に入れることができなかったので、回答として投稿してください。

于 2016-10-06T03:22:01.530 に答える
0

xとのy値をベクトルに保存し、ループの外にプロットすることをお勧めします。また、正確な解と比較するためにbigXandを出力することもできます。bigY

function y=myODE(f,y0,x0,h,x_final)
% Example:
%     f = @(x,y) (x.^2+cos(y));
%     y_final = myODE(f,0,0,0.1,2);
x=x0;
y=y0;

bigX = x0:h:x_final;
if bigX(end)<x_final
    % Doesn't occur when x_final = n*h for some integer n
    bigX = [bigX,x_final];
end
bigY = zeros(size(bigX));
count = 1;
bigY(1) = y0;

while x<x_final
  h=min(h,x_final-x);
  k1=f(x,y);
  k2=f(x+h/2,y+h*k1/2);
  k3=f(x+h/2,y+h*k2/2);
  k4=f(x+h,y+h*k3);
  y=y+(h/6)*(k1+2*k2+2*k3+k4);
  x=x+h;
  count = count+1;
  bigY(count) = y;
end

plot(bigX,bigY,'b-o')
xlabel('x')
ylabel('y')
于 2016-10-06T12:21:16.850 に答える