Manual Runge-Kutta for system of two ODEs.
37 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
Risposta accettata
Alan Stevens
il 11 Feb 2021
Modificato: Alan Stevens
il 11 Feb 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 Commenti
James Tursa
il 11 Feb 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.
Più risposte (2)
Jan
il 11 Feb 2021
Modificato: Jan
il 25 Lug 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
il 24 Lug 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 Commento
Jan
il 25 Lug 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.
Vedere anche
Categorie
Scopri di più su Numerical Integration and Differential Equations in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!