Thread: Lorentz Equations Matlab

G'day all,

I have the following Lorenz equations







I need to calculate the trajectory of in MATLAB using Heun's method (which I think is also called the trapezoidal method) with step size h=0.001, and initial conditions for t < 30.


Here is my matlab code for a single IVP (using Heun's method):


function y = heuns(f,t,x,h)
g = f(t,x);
z = x + h * g;
y = x + 0.5 * h * ( g + f(t+h,z) );


where f is a MATLAB inline function, h is the step size, t is the initial time, and x is the initial condition. However, I dont know how to use this method for the system of Lorenz equations above. Any suggestions?

I found a forum where someone called "Mr Fantastic" in a title, so I was hoping I might be able to do the same with Captainblack.

I was hoping you could help me Mr. Black, sorry, CAPTAIN Black, because you seem to have an extensive knowledge of MATLAB.

With respect to my post above, I have come up with new mfiles I think are close, but i am unsure of where I have gone wrong. Im thinking its trapezoid.m. I know its incorrect because the plot i obtain looks nothing like a lorenz equations plot should look like


function z = ydot(t,y)
% Lorenz Equations
s=10; r=28; b=8/3;
z(1) = s*y(2) - s*y(1);
z(2) = r*y(1) - y(1)*y(3) - y(2);
z(3) = y(1)*y(2)-b*y(3);

-------------------------------------------------------------------
An mfile for the trapezoid method.......(I am not 100% on this one)
-------------------------------------------------------------------
function y = trapezoid(t,x,h)
z1 = ydot(t,x);
g=x+h*z1;
z2=ydot(t+h,g);
g=x+h*z2;
z3=ydot(t+h,g);
y=x+h*(z1+z2+z3)/2;

-----------------------------------------------------------------------
function [t,y] = lorentz(int,ic,h);
n=round((int(2)-int(1))/h);
y(1, : ) =ic;
t(1) = int(1);
for k=1:n
t(k+1) = t(k) + h;
y(k+1, : ) = trapezoid(t(k),y(k,: ) ,h);
end

plot3(t,y(:,1),t,y(:,2),t,y(:,3));
-----------------------------------------------------------------
now as far as an actual matlab session,
just type:

lorentz([0 20], [5 5 5], 0.001)

Can anyone see where Ive gone wrong here?

Are you sure your algorithm is correct?

Heun's method - Wikipedia, the free encyclopedia

It seems to me it should be

Code:

function y = trapezoid(t,x,h)
yp = x + h*ydot(t,x);
y = x + h*(ydot(t,x) + ydot(t+h,yp))/2;

I tried you code, and it gives the same answer as my code. The way I wrote it is the same method, just a different way of writing it. Anyway, ive attached an image of the figure it produces.

Originally Posted by Jimmy_W

function y = trapezoid(t,x,h)
z1 = ydot(t,x);
g=x+h*z1;
z2=ydot(t+h,g);
g=x+h*z2;
z3=ydot(t+h,g);
y=x+h*(z1+z2+z3)/2;

The last line cannot be right, you have three estimates of the derivative but are dividing by 2.

RonL

Originally Posted by CaptainBlack

The last line cannot be right, you have three estimates of the derivative but are dividing by 2.

RonL

Exactly, thats why I checked it in wikipedia...It is not the same method, however it may give similar answer...

If I were you, I would check this method with linear equation.
After all these are Lorentz equations :/ we may see some "chaotic" behaviour

Point taken. Thanks, I see what i did incorrectly.

Originally Posted by albi

Exactly, thats why I checked it in wikipedia...It is not the same method, however it may give similar answer...

If I were you, I would check this method with linear equation.
After all these are Lorentz equations :/ we may see some "chaotic" behaviour

No it won't give the right answer, it is just wrong there is no way that the divisor should be 2. Just consider what it does for a constant derivative.

If this is based on a pukka method the divisor should be 3.

RonL