Can plot a function, NSolve takes too long

Question

I'm new to Mathematica, so maybe mine is an easy to solve issue, but I haven't been able to figure it out.

I have a series of linear ODEs I solve using for:

sol = ConstantArray[0, 1];
For[i = 2, i < 20, i++, 
     sol = Append[sol, 
            First[y[t] /. 
            DSolve[{y''[t] == -121/150*y'[t] - 1/10*y[t] - 
                 3 + (2/3)*D[sol[[i - 1]], t] + 1/10*sol[[i - 1]], 
                 y'[0] == 25, y''[0] == 0}, y, t]]]]

It runs fast and fine (takes a couple of seconds).

I get the first derivative of every solution:

vels = D[sol, t];

Again, everything OK... Then I plot all of them:

Plot[vels, {t, 0, 80}]

Pretty fast. I want to know the value of t for which each solution crosses a numeric boundary:

For[i = 2, i < 15, i++, 
 t10 = AppendTo[t10, 
   First[t /. NSolve[{vels[[i]] == 10, t < 40, t > 5}, t]]]]

It takes like a second or two for, say, the first 8 terms, but the 9th won't end after like an hour. It's funny because I can 'see' the result, and with a zoom large enough I could estimate the value by myself from the plot. Am I missing something important? Does it have anything to do with the methods used? Do I have to use some other command?

Thanks in advance.

Answer 1

The trouble is just that the symbolic expressions grow so large (do Last@vels and see what i mean), try this:

t10 = {}
For[i = 2, i < 15, i++, Print[i]; 
 t10 = AppendTo[t10, 
  First[t /. 
    NSolve[{Simplify[vels[[i]] // N] == 10, t < 40, t > 5}, t]]]]

(throws a few accuracy warnings..but looks right i think )

Of course even better to use functional mathematica style:

t10 = First@(t /. NSolve[{# == 10, t < 40, t > 5}, t]) & /@ 
             Simplify[vels[[2 ;; 15]] // N]

Answer 2

Note processing the solutions with N then plotting is another 'fast' approach:

Plot[N /@ vels, {t, 0, 80}]

Timing yields 0.65 seconds

Answer 3

I will use Nest instead of For, but the main suggestion I will make is to approach the problem numerically from the beginning. Even if you need the exact solution, the numerical one is fast enough that you could do both.

So first, convince yourself that the following is equivalent to your For loop:

sol2 = NestList[
   First[y[t] /. 
      DSolve[{y''[t] == -121/150*y'[t] - 1/10*y[t] - 3 + (2/3)*D[#, t] + 1/10*#,
              y'[0] == 25, y''[0] == 0}, y, t]] &,
   0, 18]; // AbsoluteTiming
(* {62.705799, Null} *)

sol == sol2
(* True *)

We can adapt this to numerical solutions with NDSolveValue (or NDSolve). It generates a warning message because y[0] is not specified and it has to calculate it (which it can do†). These are suppressed with Quiet. The output is an InterpolatingFunction instead of an expression, so I adapted how the previous solution is plugged into the current equation. Also, the initial solution must be a function, 0 &, instead of simply 0.

sol3 = Quiet[
    NestList[
     NDSolveValue[{
       y''[t] == -121/150*y'[t] - 1/10*y[t] - 3 + (2/3)*Derivative[1][#][t] + 1/10*#[t],
       y'[0] == 25, y''[0] == 0}, y, {t, 0, 80}] &,
     0 &, 18],
    NDSolveValue::icordinit]; // AbsoluteTiming
(* {0.353780, Null} *)

The output is a list of functions, so I use Through apply them to a variable t: Through[sol3[t]]. The velocities look right:

vels3 = D[Through[sol3[t]], t];

Plot[vels3, {t, 0, 80}]

Finding where the velocities are 10 takes 0.008 seconds:

(t3 = t /. FindRoot[# == 10, {t, 20, 0, 80}] & /@ Rest[vels3]) // AbsoluteTiming

(* {0.007943,
    {7.69998, 9.27022, 10.8231, 12.3619, 13.8892, 15.407, 16.9168, 18.4198, 19.9168,
     21.4087, 22.896, 24.3793, 25.8589, 27.3352, 28.8086, 30.2792, 31.7474, 33.2133}}

Comparing with the OP's vels:

(t0 = t /. FindRoot[# == 10, {t, 20, 0, 80}] & /@ Rest[vels]) // AbsoluteTiming

(* {0.533572,
    {7.69998, 9.27022, 10.8231, 12.3619, 13.8892, 15.407, 16.9168, 18.4198, 19.9168,
     21.4087, 22.896, 24.3793, 25.8589, 27.3352, 28.8086, 30.2792, 31.7474, 33.2132}} *)

Notice that the last digit of the last solution is off by 1. The accuracy of the numerical solutions can be improved with the NDSolveValue options AccuracyGoal, PrecisionGoal, and WorkingPrecision. It will take significantly longer (1 to 2 seconds, i.e. three to five times longer), but not as long as finding the exact solutions.

---

†Verification that the exact and numerical solutions agree at t -> 0.:

(sol /. t -> 0.) == (Through[sol3[t]] /. t -> 0.)
(* True *)

(There is a negligible difference in the values, a relative error of about 10^-15 at most.)