Manual Runge-Kutta for system of two ODEs.
46 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
FPixelz
el 11 de Feb. de 2021
Comentada: Jan
el 25 de Jul. de 2023
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/516862/image.jpeg)
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 comentarios
Respuesta aceptada
Alan Stevens
el 11 de Feb. de 2021
Editada: Alan Stevens
el 11 de Feb. de 2021
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 comentarios
James Tursa
el 11 de Feb. de 2021
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.
Más respuestas (2)
Jan
el 11 de Feb. de 2021
Editada: Jan
el 25 de Jul. de 2023
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);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/516882/image.png)
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.
khalida
el 24 de Jul. de 2023
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 comentario
Jan
el 25 de Jul. de 2023
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.
Ver también
Categorías
Más información sobre Numerical Integration and Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!