Solving Lorenz attractor equations using Runge kutta (RK4) method
59 views (last 30 days)
Show older comments
reema shrestha
on 11 Oct 2017
Commented: Gerardo ramirez
on 18 May 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 Comments
Accepted Answer
Mischa Kim
on 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 Comments
More Answers (1)
tyfcwgl
on 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)
0 Comments
See Also
Categories
Find more on Digital Filter Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!