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

1 Answer 1

Reset to default
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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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