1

I am trying to simulate a coupled system of N blocks and springs. I am using NDSolve. For N = 25 everything goes smoothly, however for N > 25 it doesn't work. No error message, no warning, just the usual beep.

Is there any limitation on the number of ODEs NDSolve can handle? (I'm using V9.0.1.0, Home Edition.)

UPDATE: The problem is not related with some kind of restriction on NDSolve, I just checked it with a more simpler system and n=500; NDSolve worked fine. What follows is the code; as I said before for n=25 it works well, in less than a second, however for n>26 It doesn't. Why the beep? says : The kernel Local has quit (exited) during the course of an evaluation

n = 26;
U0 = 20;
Nt = 3; 
b = 0.5;
umax = 2;
tr = 40.0;


S4 = Join[{u[1]''[t] == 
 1/Nt (U0*(1 - Exp[-(t/tr)]) - u[1][t]) - 
  2*u[1][t]*
   HeavisideTheta[2*umax - Abs[2*u[1][t]]] + (u[n][t] - 
    u[1][t]) + (u[2][t] - u[1][t]) - b*u[1]'[t], 
u[2]''[t] == 
 1/Nt (U0*(1 - Exp[-(t/tr)]) - u[2][t]) - 
  2*u[2][t]*
   HeavisideTheta[2*umax - Abs[2*u[2][t]]] + (u[2 - 1][t] - 
    u[2][t]) + (u[2 + 1][t] - u[2][t]) - b*u[2]'[t],
u[3]''[t] == 
 1/Nt (U0*(1 - Exp[-(t/tr)]) - u[3][t]) + (u[3 - 1][t] - 
    u[3][t]) + (u[3 + 1][t] - u[3][t]) - b*u[3]'[t]
}, Table[
u[i]''[t] == 
 1/Nt (U0*(1 - Exp[-(t/tr)]) - u[i][t]) - 
  2*u[i][t]*
   HeavisideTheta[2*umax - Abs[2*u[i][t]]] + (u[i - 1][t] - 
    u[i][t]) + (u[i + 1][t] - u[i][t]) - b*u[i]'[t], {i, 4, 
 n - 1}],
{u[n]''[t] == 
 1/Nt (U0*(1 - Exp[-(t/tr)]) - u[n][t]) - 
  2*u[n][t]*
   HeavisideTheta[2*umax - Abs[2*u[n][t]]] + (u[1][t] - 
    u[n][t]) + (u[n - 1][t] - u[n][t]) - b*u[n]'[t]},
Table[u[i][0] == 0, {i, 1, n}],
Table[u[i]'[0] == 0, {i, 1, n}]];
Sol = NDSolve[S4, Table[u[i], {i, 1, n}], {t, 0, 70}, 
Method -> "ExplicitRungeKutta"];
Plot[Evaluate[Table[u[i][t], {i, 1, n}] /. Sol], {t, 0, 70}]
Jorge
  • 51
  • 6
  • 3
    You may be running out of memory. – bbgodfrey Dec 31 '15 at 10:03
  • Please provide a minimal example that can reproduce the issue. – xzczd Dec 31 '15 at 10:14
  • For example: 
    n=10;$$ $$

func = ToExpression@Table["f" <> ToString[i] <> "[x]", {i, 0, n}]; $$ $$


NDSolve[Flatten@{D[#[[1]], x] == #[[2]] & /@ 
 Transpose@{func, RotateLeft@func}, (# /. x -> 0) == 1 & /@ func},

    func, {x, 0, 1}] // Short $$ $$

    works fine for any number of ODE's (I checked up to $n=1000$). This system of ODE's is $f'_{i-1}(x)=f_i(x)$, with $f_i(0)=1$.

     – AccidentalFourierTransform Dec 31 '15 at 16:57
  • Have you tried the menu command, Help -> Why the Beep?...? – Michael E2 Jan 01 '16 at 01:36
  • @AccidentalFourierTransform You're right, I just checked it with a more simpler coupled system of ODEs and NDSolve worked fine for N=500. So the problem is not the number of equations. – Jorge Jan 01 '16 at 07:49
  • @MichaelE2 Why the beep? shows: The kernel Local has quit (exited) during the course of an evaluation I'm gonna check that. – Jorge Jan 01 '16 at 07:57
  • Instead of writing out the same equation multiple times, write it in vector form. http://mathematica.stackexchange.com/questions/78641/vector-form-using-ndsolve – Szabolcs Jan 01 '16 at 10:33

1 Answers1

3

Finally I realized that the problem is not related directly with NDSolve but with RungeKutta method and HeavisideTheta function. I changed the latter by a more "softer" function:

(1 + Tanh[p (1 - (x/xmax))])/(1 + Tanh[p])

where p is a number that measures the stiffness of the function and xmax is approximately the zone where the function goes to zero. Now I can work with hundreds of ODE´s with no problem.

Jorge
  • 51
  • 6