MATLAB как мне построить мое решение ОДУ Рунге-Kutta4 на всем интервале

Question

я сделал этот код для решений ОДУ с использованием Рунге-Kutta4:

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 является функция у»= F (х, у), y0 - начальное значение, x0 - это функция, начинающаяся с h и x_final, где функция останавливается.

Я попробовал свой код и он решает ОДУ для меня правильно, но я также хочу, чтобы построить свою функцию по х оси на интервале x0 к x_final с h подынтервалами. Когда я пытаюсь построить его, используя plot(x0:h:x_final,y), я получаю пустой граф. Я понимаю (догадываясь), что мне нужно привязать мой y к нескольким x, чтобы построить, но как я могу это сделать, не меняя слишком много кода?

Как я могу построить график для y заданного y0, интервал x0 к x_final и дал h?

Новое в MATLAB, оцените всю помощь, которую я могу получить!

Редактировать: Чтобы уточнить, для чего предназначен мой код;

Мне нужен этот ODE-решатель как для решения, так и для графического отображения. Я должен изучить ошибку усечения, посмотрев на значения y на h по сравнению с 2*h и стабильность Runge-Kutta4, посмотрев на графики y с различными h.

Answer 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 

Может быть, завтра я буду переписывать этот ответ на что-то правильное, но — для теперь — доброй ночи! :-)

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')