如何在整个区间上绘制龙格-库塔4 ODE解

Question

我使用runge - kutt4为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到x_final的区间上用h子区间在xy轴上绘制我的函数。当我尝试使用plot(x0:h:x_final,y)绘制它时，我只得到一个空图。我理解(猜测)，我必须绑定我的y到几个x为了绘图，但我怎么能做到这一点，而不改变我的代码太多?

我如何绘制y给定y0，区间x0到x_final和给定h的图形?

新的MATLAB，感谢所有的帮助，我可以得到!

编辑:明确我的代码是什么;

我需要这个ODE求解器来求解和绘图。我应该通过查看h上y与2*h的值来研究截断误差，并通过查看y与不同h的图形来研究Runge-Kutta4的稳定性。

Answer 1

这不是一个非常聪明的代码重构(因为它会减慢求解速度，还会破坏你的图形，这取决于你在ODE中有多少步骤)，但我很困，所以我选择了hack:

    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

也许明天我会把这个答案重新写一遍，但是现在，晚安!: -)

Answer 2

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

评论框不允许我将固定代码放在"代码格式"中，所以将其作为答案发布。

Answer 3

我建议您将x和y值存储到向量中，并在循环外绘制它们。您可能还想输出bigX和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')