Manual Runge-Kutta for system of two ODEs.
37 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponse acceptée
Alan Stevens
le 11 Fév 2021
Modifié(e) : Alan Stevens
le 11 Fév 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 commentaires
James Tursa
le 11 Fév 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.
Plus de réponses (2)
Jan
le 11 Fév 2021
Modifié(e) : Jan
le 25 Juil 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
le 24 Juil 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 commentaire
Jan
le 25 Juil 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.
Voir également
Catégories
En savoir plus sur Numerical Integration and Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!