3

I have defined a function which takes a coefficient for a differential equation as an argument, solves the differential equation with NDSolve, and returns the value of the solution at time t=10. I can plot this, but I can't seem to find the minimum value.

This code demonstrates the problem i'm having:

f[k_] := y[10] /. NDSolve[
    {y''[t] == -k*y[t], y[1] == 2 , y'[2] == 1},
    y,
    {t, 0, 10}
];
Plot[f[k], {k, 0, 2}] (* works great *)
NMinimize[{f[k], k > 0, k < .75}, k] (* fails *)

The NMinimize fails with the error NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 1.'.

So how do I find the minimum of my function f?

CesiumLifeJacket
  • 33
  • 3
  • 2
    Try f[k_?NumericQ] := y[10] /. First @ NDSolve[{y''[t] == -k*y[t], y[1] == 2, y'[2] == 1}, y, {t, 0, 10}]. See the pitfalls question. You also might want to consider FindMinimum instead of NMinimize. – Michael E2 May 08 '14 at 23:36
  • Be sure to Clear[f] before you try Michael E2's solution, otherwise Mathematica will overload f with two definitions and use the less-restrictive one that you don't want. – Chris K May 08 '14 at 23:38

2 Answers2

4

Although Michael E2's suggestion in his comment about using the argument pattern k_?NumericQ will certainly work, I would like to suggest a double interpolation approach because it is much faster.

ff[k_] := NDSolve[{y''[t] == -k*y[t], y[1.] == 2, y'[2.] == 1.}, y, {t, 0., 10.}][[1, 1, 2]]
AbsoluteTiming[
  f = Interpolation@Table[{k, ff[k][10.]}, {k, 0., 2, .01}];]

{2.748679, Null}

AbsoluteTiming @ NMinimize[{f[k], {k > 0., k < .5}}, k]

{0.075349, {-3.95571, {k -> 0.212177}}}

Now compare this with 

g[k_?NumericQ] := 
  NDSolve[{y''[t] == -k*y[t], y[1.] == 2, y'[2.] == 1.}, y, {t, 0., 10.}][[1, 1, 2]][10.]
AbsoluteTiming @ NMinimize[{g[k], {k > 0., k < .5}}, k]

{44.634075, {-3.95571, {k -> 0.212174}}}

The first result for the minimum is a little less precise, but much faster.

m_goldberg
  • 107,779
  • 16
  • 103
  • 257
  • I think FindMinimum is even faster. On this particular problem, it's even faster with DSolve. – Michael E2 May 09 '14 at 01:02
  • @MichaelE2. Didn't get that far out the box when I was thinking about out-of-box ideas to improve the wretched performance of the OP approach :) – m_goldberg May 09 '14 at 02:03
  • @m_goldberg could you explain a little why your double interpolation approach works out so fast, and how broadly it can be applied? – Chris K May 09 '14 at 02:22
  • 1
    @ChrisK. I would say it's because ff in my method gets called many fewer times than g does in my version of the OP's method. But in this case, you should look at J. W. Perry's answer which shows an even better approach. However, it is a generally applicable trick and you might find other occasions where it will help. – m_goldberg May 09 '14 at 02:50
2

For the record, here is my implementation of DSolve and FindMinimum on the same equation. This does not address your NDSolve work but the result is relatively instantaneous.

soln = DSolve[{y''[t] == -k*y[t], y[1] == 2, y'[2] == 1}, y[t], t];
f[k_] := Evaluate[y[t] /. soln] /. t -> 10
FindMinimum[{f[k][[1]], 0 <= k <= 3/4}, {k, 1/10}]

I plot the same with

Plot[f[k], {k, 0, 2}]
J. W. Perry
  • 1,080
  • 1
  • 6
  • 11