Solving Lorenz attractor equations using Runge kutta (RK4) method
63 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
reema shrestha
am 11 Okt. 2017
Kommentiert: Gerardo ramirez
am 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 Kommentare
Akzeptierte Antwort
Mischa Kim
am 11 Okt. 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 Kommentare
Weitere Antworten (1)
tyfcwgl
am 29 Okt. 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 Kommentare
Siehe auch
Kategorien
Mehr zu Digital Filter Analysis finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!