# Problems with solving two coupled lorenz equation with NDsolve

I have a system of two coupled lorenz equations and I'm to integrate it:

coupledlorenz =
NDSolve[{Subscript[x, 1]'[t] ==
10*(Subscript[x, 1][t] - Subscript[x, 2][t]),
Subscript[x, 2]'[t] ==
40 Subscript[x, 1][t] - Subscript[x, 2][t] -
Subscript[x, 1][t] Subscript[x, 3][t],
Subscript[x, 3]'[t] ==
Subscript[x, 1][t] Subscript[x, 2][t] - (8/3) Subscript[x, 3][t],
Subscript[y, 1]'[t] ==
10 (Subscript[y, 2][t] - Subscript[y, 1][t]) +
0.1 (Subscript[x, 1][t] - Subscript[y, 1][t]),
Subscript[y, 2]'[t] ==
35 Subscript[y, 1][t] - Subscript[y, 2][t] -
Subscript[y, 1][t] Subscript[y, 2][t],
Subscript[y, 3]'[t] ==
Subscript[y, 1][t] Subscript[y, 2][t] - (8/3) Subscript[y, 3][t],
Subscript[x, 1][0] == 0, Subscript[x, 2][0] == 1,
Subscript[x, 3][0] == 0, Subscript[y, 1][0] == 0,
Subscript[y, 2][0] == 1, Subscript[y, 3][0] == 0}, {Subscript[x,
1], Subscript[x, 2]}, {t, 0, 500}, MaxSteps -> Infinity]


The problem is that Mathematica (Version 10.1) is running and running and eats up a lot of memory. I presume the problem is with the MaxSteps-> Infinity option since reducing MaxSteps works. On the other hand in the documentation it is recommended to use MaxSteps->Infinity for a lorenz system.

In short, those equations are from a book where they solved this system with a fourth-order Runge-Kutta with time step 0.03 and t=10000. So it tried it like this after browsing through some of the other threads:

    ClassicalRungeKuttaCoefficients[4, prec_] :=
With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}},
bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}},
N[{amat, bvec, cvec}, prec]]

coupledlorenz =
NDSolve[{Subscript[x, 1]'[t] ==
10*(Subscript[x, 1][t] - Subscript[x, 2][t]),
Subscript[x, 2]'[t] ==
40 Subscript[x, 1][t] - Subscript[x, 2][t] -
Subscript[x, 1][t] Subscript[x, 3][t],
Subscript[x, 3]'[t] ==
Subscript[x, 1][t] Subscript[x, 2][t] - (8/3) Subscript[x, 3][t],
Subscript[y, 1]'[t] ==
10 (Subscript[y, 2][t] - Subscript[y, 1][t]) +
1 (Subscript[x, 1][t] - Subscript[y, 1][t]),
Subscript[y, 2]'[t] ==
35 Subscript[y, 1][t] - Subscript[y, 2][t] -
Subscript[y, 1][t] Subscript[y, 2][t],
Subscript[y, 3]'[t] ==
Subscript[y, 1][t] Subscript[y, 2][t] - (8/3) Subscript[y, 3][t],
Subscript[x, 1][0] == 0, Subscript[x, 2][0] == 1,
Subscript[x, 3][0] == 0, Subscript[y, 1][0] == 0,
Subscript[y, 2][0] == 1, Subscript[y, 3][0] == 0}, {Subscript[x,
1], Subscript[x, 2], Subscript[x, 3]}, {t, 0, 500},
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,
"Coefficients" -> ClassicalRungeKuttaCoefficients},
StartingStepSize -> 0.03]


This doesn't work either and gives me the following error:

General::ovfl: Overflow occurred in computation. >>

General::ovfl: Overflow occurred in computation. >>

NDSolve::nlnum: The function value {-4.153140717608156*10^1316669612242674,Overflow[],<<3>>,-2.511471835313141*10^1316669612242539} is not a list of numbers with dimensions {6} at {t,Subscript[x, 1][t],Subscript[x, 2][t],Subscript[x, 3][t],Subscript[y, 1][t],Subscript[y, 2][t],Subscript[y, 3][t]} = {1.755,3.848459941679012*10^813746572320117,<<4>>,7.073363386337455*10^813746572319983}. >>


I just don't understand why it is not working and how I can integrate this system. It must be integrable, since the authors of the book used the same system, so I must be doing something totally wrong here.

Any help appreciated.

Never use Subscript when debugging Mathematica code, because it is difficult to read. The code can be rewritten for clarity as,

coupledlorenz = NDSolve[{
x1'[t] == 10*(x1[t] - x2[t]),
x2'[t] == 40 x1[t] - x2[t] - x1[t] x3[t],
x3'[t] == x1[t] x2[t] - (8/3) x3[t],
y1'[t] == 10 (y2[t] - y1[t]) + 0.1 (x1[t] - y1[t]),
y2'[t] == 35 y1[t] - y2[t] - y1[t] y2[t],
y3'[t] == y1[t] y2[t] - (8/3) y3[t],
x1[0] == 0, x2[0] == 1, x3[0] == 0, y1[0] == 0, y2[0] == 1,
y3[0] == 0},
{x1, x2, x3, y1, y2, y3}, {t, 0, 2}];


which runs to completion in a minute or so. Values at t = 2 are

Flatten[{x1[2], x2[2], x3[2], y1[2], y2[2], y3[2]} /. coupledlorenz]
(* {646430., 4.76559, 38.6208, 3250.67, 34.9892, 9315.7} *)


Clearly, x1 is growing rapidly. To gain more insight, note that the x variables do not depend on the y variables, so the first three equations can be solved separately for simplicity.

s = NDSolve[{
x1'[t] == 10*(x1[t] - x2[t]),
x2'[t] == 40 x1[t] - x2[t] - x1[t] x3[t],
x3'[t] == x1[t] x2[t] - (8/3) x3[t],
x1[0] == 0, x2[0] == 1, x3[0] == 0},
{x1, x2, x3}, {t, 0, 2}];
Plot[Evaluate[{x1[t], x2[t], x3[t]} /. s], {t, 0, 2}]


x1 is growing exponentially, and examination of the equations indicates that its value is proportional to Exp[10 t] for large t. x2 on the other hand is gradually coverging to 0, and x3 to 40. Moreover, their oscillation rates are given roughly by x1. In other words, they are oscillating enormously fast as t increases, and for this reason the time step must decrease as 1/x1. Not surprisingly integrating to large t takes seemingly forever.

• Thank you very much, that makes it much more clearer for me. Won't use subscripts again, thanks for the hint :)! But I still don't understand how the authors of the book where i found the equations could use 10 000 time points for the integration of this system. Maybe I understood somethings wrong :(..anyway.. thanks again :) – holistic Aug 15 '15 at 15:29