Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Mathematica Understanding this Computer program

  1. Feb 26, 2015 #1
    Below is a computer program for Mathematica that details how to calculate the Lyapunov exponent of a differential equation. What is the purpose of the code colored in red?

    ClearAll["Global`*"];
    deq1 =-(y1[t]+ z1[t]);
    deq2 = x1[t]+0.1 y1[t];
    deq3 =0.2+ x1[t] z1[t]-5.7 z1[t];

    deq4 =-(y2[t]+ z2[t]);
    deq5 = x2[t]+0.1 y2[t];
    deq6 =0.2+ x2[t] z2[t]-5.7 z2[t];

    x10 =1; y10 =1; z10 =1;
    dx0 =10^-8;
    x20 = x10 + dx0; y20 = y10; z20 = z10;
    tin =0; tfin =10000;
    tstep =1;
    acc =12;

    lcedata ={};
    sum =0;

    d0 =Sqrt[(x10 - x20)^2+(y10 - y20)^2+(z10 - z20)^2];

    For[i =1, i < tfin/tstep, i++,

    sdeq ={x1'[t]== deq1, y1'[t]== deq2, z1'[t]== deq3,
    x2'[t]== deq4, y2'[t]== deq5, z2'[t]== deq6, x1[0]== x10,
    y1[0]== y10, z1[0]== z10, x2[0]== x20, y2[0]== y20,
    z2[0]== z20};
    sol =NDSolve[
    sdeq,{x1[t], y1[t], z1[t], x2[t], y2[t], z2[t]},{t,0, tstep},MaxSteps->Infinity,Method->"Adams",PrecisionGoal-> acc,AccuracyGoal-> acc];

    xx1[t_]= x1[t]/. sol[[1]];
    yy1[t_]= y1[t]/. sol[[1]];
    zz1[t_]= z1[t]/. sol[[1]];

    xx2[t_]= x2[t]/. sol[[1]];
    yy2[t_]= y2[t]/. sol[[1]];
    zz2[t_]= z2[t]/. sol[[1]];

    d1 =Sqrt[(xx1[tstep]- xx2[tstep])^2+(yy1[tstep]- yy2[tstep])^2+(zz1[tstep]- zz2[tstep])^2];

    sum +=Log[d1/d0];
    dlce = sum/(tstep*i);AppendTo[lcedata,{tstep*i,Log10[dlce]}];

    w1 =(xx1[tstep]- xx2[tstep])*(d0/d1);
    w2 =(yy1[tstep]- yy2[tstep])*(d0/d1);
    w3 =(zz1[tstep]- zz2[tstep])*(d0/d1);

    x10 = xx1[tstep];
    y10 = yy1[tstep];
    z10 = zz1[tstep];

    x20 = x10 + w1;
    y20 = y10 + w2;
    z20 = z10 + w3;

    i = i++;


    If[Mod[tstep*i,100]==0,Print[" For t = ", tstep*i," , "," LCE = ", dlce]]

    ]

    S0 =ListPlot[{lcedata},Frame->True,Axes->False,PlotRange->All,Joined->True,FrameLabel->{"t","log10(LCE)"},FrameStyle->Directive["Helvetica",17],ImageSize->550]




    The link to the code: http://mathematica.stackexchange.com/questions/17593/lyapunov-exponent/17608#17608
     
  2. jcsd
  3. Feb 26, 2015 #2

    jedishrfu

    Staff: Mentor

    This is how you numerically integrate a function computationally.

    As an example, say you want to integrate dy=2xdx over some interval of 0 to 5

    Then you'd do something like this:

    Code (Text):

    dx = 0.01
    x = 0
    y = 0
    For i from 0 to 500
      x = x + dx
      dy = 2 * x * dx
      y = y + dy
      Print x,y
    Next
     
     
  4. Feb 26, 2015 #3
    I don't think that answers my question.
     
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook




Loading...