1
$\begingroup$

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}]
$\endgroup$
  • 3
    $\begingroup$ You may be running out of memory. $\endgroup$ – bbgodfrey Dec 31 '15 at 10:03
  • $\begingroup$ Please provide a minimal example that can reproduce the issue. $\endgroup$ – xzczd Dec 31 '15 at 10:14
  • $\begingroup$ 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$. $\endgroup$ – AccidentalFourierTransform Dec 31 '15 at 16:57
  • $\begingroup$ Have you tried the menu command, Help -> Why the Beep?...? $\endgroup$ – Michael E2 Jan 1 '16 at 1:36
  • $\begingroup$ @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. $\endgroup$ – Jorge Jan 1 '16 at 7:49
3
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.