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

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?

    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


    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
  4. Feb 26, 2015 #3
    I don't think that answers my question.
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Similar Discussions: Understanding this Computer program
  1. Programming in Matlab (Replies: 4)