Results 1 to 8 of 8

Math Help - Lorentz Equations Matlab

  1. #1
    Junior Member
    Joined
    Sep 2008
    Posts
    47

    Lorentz Equations Matlab

    G'day all,

    I have the following Lorenz equations

     x' = 10(y - x)

     y' = x(28 - z) -y

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

    I need to calculate the trajectory of  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  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?
    Follow Math Help Forum on Facebook and Google+

  2. #2
    Junior Member
    Joined
    Sep 2008
    Posts
    47

    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.

    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?
    Last edited by Jimmy_W; September 25th 2008 at 09:17 AM. Reason: forgot something
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Member
    Joined
    Jan 2006
    From
    Gdansk, Poland
    Posts
    117
    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;
    Follow Math Help Forum on Facebook and Google+

  4. #4
    Junior Member
    Joined
    Sep 2008
    Posts
    47
    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.
    Attached Thumbnails Attached Thumbnails Lorentz Equations Matlab-lorenz-using-trap.jpg  
    Follow Math Help Forum on Facebook and Google+

  5. #5
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Quote Originally Posted by Jimmy_W View Post
    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
    Follow Math Help Forum on Facebook and Google+

  6. #6
    Member
    Joined
    Jan 2006
    From
    Gdansk, Poland
    Posts
    117
    Quote Originally Posted by CaptainBlack View Post
    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
    Follow Math Help Forum on Facebook and Google+

  7. #7
    Junior Member
    Joined
    Sep 2008
    Posts
    47
    Point taken. Thanks, I see what i did incorrectly.
    Follow Math Help Forum on Facebook and Google+

  8. #8
    Grand Panjandrum
    Joined
    Nov 2005
    From
    someplace
    Posts
    14,972
    Thanks
    4
    Quote Originally Posted by albi View Post
    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
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. Equilibrium Points in the Lorentz System
    Posted in the Algebra Forum
    Replies: 1
    Last Post: August 11th 2013, 11:18 PM
  2. Special Lorentz Transformation Question
    Posted in the Advanced Math Topics Forum
    Replies: 1
    Last Post: November 14th 2010, 10:10 AM
  3. lorentz transformation question..
    Posted in the Advanced Applied Math Forum
    Replies: 2
    Last Post: January 1st 2010, 01:50 AM
  4. Lorentz force law
    Posted in the Differential Equations Forum
    Replies: 1
    Last Post: September 20th 2009, 01:54 AM
  5. lorentz transformations
    Posted in the Advanced Applied Math Forum
    Replies: 0
    Last Post: January 23rd 2008, 03:11 PM

Search Tags


/mathhelpforum @mathhelpforum