Solving Lorenz attractor equations using Runge kutta (RK4) method
21 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
reema shrestha
le 11 Oct 2017
Commenté : Gerardo ramirez
le 18 Mai 2020
I am trying to write a code for the simulation of lorenz attractor using rk4 method. Here is the code:
clc;
clear all;
t(1)=0; %initializing x,y,z,t
x(1)=1;
y(1)=1;
z(1)=1;
sigma=10; %value of constants
rho=28;
beta=(8/3);
h=0.1; %step size
t=0:h:100;
f=@(t,x,y,z) sigma*(y-x); %ode
g=@(t,x,y,z) x*rho-x.*z-y;
p=@(t,x,y,z) x.*y-beta*z;
for i=1:(length(t)-1) %loop
k1=h*f(t(i),x(i),y(i),z(i));
l1=h*g(t(i),x(i),y(i),z(i));
m1=h*p(t(i),x(i),y(i),z(i));
k2=h*f(t(i)+h,(x(i)+0.5*k1),(y(i)+(0.5*l1)),(z(i)+(0.5*m1)));
l2=h*f(t(i)+h,(x(i)+0.5*k1),(y(i)+(0.5*l1)),(z(i)+(0.5*m1)));
m2=h*f(t(i)+h,(x(i)+0.5*k1),(y(i)+(0.5*l1)),(z(i)+(0.5*m1)));
k3=h*f(t(i)+h,(x(i)+0.5*k2),(y(i)+(0.5*l2)),(z(i)+(0.5*m2)));
l3=h*f(t(i)+h,(x(i)+0.5*k2),(y(i)+(0.5*l2)),(z(i)+(0.5*m2)));
m3=h*f(t(i)+h,(x(i)+0.5*k2),(y(i)+(0.5*l2)),(z(i)+(0.5*m2)));
k4=h*f(t(i)+h,(x(i)+k3),(y(i)+l3),(z(i)+m3));
l4=h*g(t(i)+h,(x(i)+k3),(y(i)+l3),(z(i)+m3));
m4=h*p(t(i)+h,(x(i)+k3),(y(i)+l3),(z(i)+m3));
x(i+1)=x(i)+h*(1/6)*(k1+2*k2+2*k3+k4); %final equations
y(i+1)=y(i)+h*(1/6)*(k1+2*k2+2*k3+k4);
z(i+1)=z(i)+h*(1/6)*(m1+2*m2+2*m3+m4);
end
plot3(x,y,z)
But the solutions are not right. I don't know what to do. I know we can do using ode solvers but i wanted to do using rk4 method. I searched for the solutions in different sites but i didn't find many using rk4. While there were some but only algorithm. I tried to compare my solutions with ode45 but doesn't match at all. it's totally different.
0 commentaires
Réponse acceptée
Mischa Kim
le 11 Oct 2017
The only bug that I can see at first glance is here
y(i+1) = y(i) + h*(1/6)*(l1+2*l2+2*l3+l4); % replace ki by li
You also might want to play with (decrease) the step size.
5 commentaires
Plus de réponses (1)
tyfcwgl
le 29 Oct 2017
I think there are many bugs in your code. After my modifying, it works well. the code and result are below.
clc;
clear all;
t(1)=0; %initializing x,y,z,t
x(1)=1;
y(1)=1;
z(1)=1;
sigma=10; %value of constants
rho=28;
beta=(8.0/3.0);
h=0.01; %step size
t=0:h:20;
f=@(t,x,y,z) sigma*(y-x); %ode
g=@(t,x,y,z) x*rho-x.*z-y;
p=@(t,x,y,z) x.*y-beta*z;
for i=1:(length(t)-1) %loop
k1=f(t(i),x(i),y(i),z(i));
l1=g(t(i),x(i),y(i),z(i));
m1=p(t(i),x(i),y(i),z(i));
k2=f(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));
l2=g(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));
m2=p(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));
k3=f(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
l3=g(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
m3=p(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
k4=f(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
l4=g(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
m4=p(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
x(i+1) = x(i) + h*(k1 +2*k2 +2*k3 +k4)/6; %final equations
y(i+1) = y(i) + h*(l1 +2*l2 +2*l3 +l4)/6;
z(i+1) = z(i) + h*(m1+2*m2 +2*m3 +m4)/6;
end
plot3(x,y,z)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/178987/image.png)
0 commentaires
Voir également
Catégories
En savoir plus sur Digital Filter Analysis 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!