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，但我也想在 xy 轴上绘制我的函数到子间隔x0的x_final区间h。当我尝试使用它来绘制它时，plot(x0:h:x_final,y)我只得到一个空图。我知道（猜测）我必须将我的绑定y到几个x才能进行绘图，但是我怎样才能在不过多更改代码的情况下做到这一点？

如何绘制ygiven y0、 interval x0tox_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

我建议您将xandy值存储到向量中并将它们绘制在循环之外。您可能还想输出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')
于 2016-10-06T12:21:16.850 回答
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 回答