Manual Runge-Kutta for system of two ODEs.

Question

FPixelz 2021-2-11

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/742732-manual-runge-kutta-for-system-of-two-odes

评论： Jan 2023-7-25

采纳的回答： Alan Stevens

在 MATLAB Online 中打开

I am struggling to obtain the correct graph for the system of ODEs as follows:

x'=-y+6x, y'=-y+4x, between t=0,0.7

I can obtain the correct graph using Euler's method, as seen here:

But cannot do the same for a manual Runge Kutta method. And I don't want to use the integrated ode45 functions if I don't have to. What am I doing wrong? My code is below:

clear,clc
h = 0.1
t_beg = 0
t_end = 0.7
x_initial= 0.5
y_initial= -0.5
F_tx=@(x,y)(-x+6*y);
F_ty=@(x,y)(-y+4*x);
t=t_beg:h:t_end;
x=zeros(1,length(t));
x(1)=x_initial;
y=zeros(1,length(t));
y(1)=y_initial;
for i=1:(length(t)-1)
    
    kx1 = F_tx(t(i),x(i));
    kx2 = F_tx(t(i)+0.5*h,x(i)+0.5*h*kx1);
    kx3 = F_tx((t(i)+0.5*h),(x(i)+0.5*h*kx2));
    kx4 = F_tx((t(i)+h),(x(i)+kx3*h));
    
    x(i+1) = x(i) + (1/6)*(kx1+2*kx2+2*kx3+kx4)*h;
    
    ky1 = F_ty(t(i),y(i));
    ky2 = F_ty(t(i)+0.5*h,y(i)+0.5*h*ky1);
    ky3 = F_ty((t(i)+0.5*h),(y(i)+0.5*h*ky2));
    ky4 = F_ty((t(i)+h),(y(i)+ky3*h));
    
    y(i+1) = y(i) + (1/6)*(ky1+2*ky2+2*ky3+ky4)*h;
    
end
% plot(x,y)
figure(1)
plot(t,y)
hold on
plot(t,x)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2021-2-11

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/742732-manual-runge-kutta-for-system-of-two-odes#answer_621117

编辑：Alan Stevens 2021-2-11

在 MATLAB Online 中打开

You need to change the order within the loop to

for i=1:(length(t)-1)
    
    kx1 = F_tx(x(i),y(i));
    ky1 = F_ty(x(i),y(i));
    kx2 = F_tx(x(i)+0.5*h*kx1,y(i)+0.5*h*ky1);
    ky2 = F_ty(x(i)+0.5*h*kx1,y(i)+0.5*h*ky1);
    kx3 = F_tx((x(i)+0.5*h*kx2),y(i)+0.5*h*ky2);
    ky3 = F_ty((x(i)+0.5*h*kx2),(y(i)+0.5*h*ky2));
    kx4 = F_tx((x(i)+kx3*h),y(i)+ky3*h);
    ky4 = F_ty((x(i)+kx3*h),(y(i)+ky3*h));
    
    x(i+1) = x(i) + (1/6)*(kx1+2*kx2+2*kx3+kx4)*h; 
    y(i+1) = y(i) + (1/6)*(ky1+2*ky2+2*ky3+ky4)*h;
    
end

and note that the first argument is x not t.

4 个评论
显示 2更早的评论隐藏 2更早的评论

James Tursa 2021-2-11

For three variables x, y, z you still need to respect the order of the k evaluations. Do kx1, ky1, kz1 first. Then do kx2, ky2, kz2. Etc.

Or to get this same effect use the vector approach that Jan has posted.

FPixelz 2021-2-16

Thanks for your help!

请先登录，再进行评论。

Answer 2

Jan 2021-2-11

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/742732-manual-runge-kutta-for-system-of-two-odes#answer_621127

编辑：Jan 2023-7-25

在 MATLAB Online 中打开

The diagram looks, yoike you are integrating:

@(t, y) [-y(1)+6*y(2); -y(2)+4*y(1)]

This code produces an equivalent output:

t0 = 0;
tF = 0.7;
x0 = 0.5;
y0 = -0.5;
[t, y] = ode45(@(t,y) [-y(1)+6*y(2); -y(2)+4*y(1)], ...
               [t0, tF], [x0, y0]);
figure;
plot(t,y);

But your function to be integrated is something else:

F_tx = @(x,y) (-x + 6 * y);
F_ty = @(x,y) (-y + 4 * x);

Here the function depends on the 1st input, which is t in my code. This is a confusion of "x/y" versus "t/y", whereby your "y" consists of the components x and y.

[EDITED] A working solution:

F = @(t, y) [-y(1) + 6 * y(2); ...
             -y(2) + 4 * y(1)];
y       = zeros(2, length(t));
y(:, 1) = [x_initial; y_initial];
for i=1:(length(t)-1)    
    kx1 = F(t(i),           y(:, i));
    kx2 = F(t(i) + 0.5 * h, y(:, i) + 0.5 * h * kx1);
    kx3 = F(t(i) + 0.5 * h, y(:, i) + 0.5 * h * kx2);
    kx4 = F(t(i) + h,       y(:, i) + kx3 * h);
    
    y(:, i+1) = y(:, i) + (kx1 + 2 * kx2 + 2 * kx3 + kx4) * h / 6;
end
figure()
plot(t, y)

By the way, compare the readability of the code, which contains spaces around the operators.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

FPixelz 2021-2-16

Thanks for the help, this has proved very useful to me.

请先登录，再进行评论。

Answer 3

khalida 2023-7-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/742732-manual-runge-kutta-for-system-of-two-odes#answer_1277723

I am struggling to obtain the correct algorithm for the system of ODEs as follows:

y'=z, z'=-((1+5x)/(2x(x+1))z-y^3+5x+11x^2+0.296x^9+0.666x^10+0.5x^11+0.125x^12), between x=0,0.7

y(0)=0, y'(0)=z(0)=0

i need the matlab code for this system of odes

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jan 2023-7-25

Please do not append a new question in the section for answers of another question. Open a new thread instead and delete this message here. Post, what you have tried so far.

请先登录，再进行评论。