Lorentz Equations Matlab

• Sep 24th 2008, 06:01 PM
Jimmy_W
Lorentz Equations Matlab
G'day all,

I have the following Lorenz equations

$\displaystyle x' = 10(y - x)$

$\displaystyle y' = x(28 - z) -y$

$\displaystyle z' = xy - \frac{8}{3}z$

I need to calculate the trajectory of $\displaystyle Y(t) = (x(t), y(t), z(t))$ in MATLAB using Heun's method (which I think is also called the trapezoidal method) with step size h=0.001, and initial conditions $\displaystyle Y(0) = (5,5,5)$ 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?
(Worried)
• Sep 25th 2008, 08:54 AM
Jimmy_W
Calling CaptainBlack
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. (Itwasntme)

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(Worried)

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?
• Sep 25th 2008, 02:37 PM
albi
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;
• Sep 25th 2008, 07:27 PM
Jimmy_W
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.
• Sep 25th 2008, 10:43 PM
CaptainBlack
Quote:

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
• Sep 26th 2008, 03:20 AM
albi
Quote:

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
• Sep 26th 2008, 03:31 AM
Jimmy_W
Point taken. Thanks, I see what i did incorrectly.
• Sep 26th 2008, 04:04 AM
CaptainBlack
Quote:

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