Maximising solutions of NDSolve

Question

I have the following system of differential equations

systeqns[gn_, gp_, sn_, sp_, K_, a_, 
  b_, n0_, p0_, tg_] := {

  n'[t] == 
   If[t <= tg, gn, sn]*n[t]*(1 - (n[t] + p[t])/K) - 
    a*n[t] + b*p[t],

  p'[t] == 
   If[t <= tg, gp, sp]*p[t]*(1 - (n[t] + p[t])/K) + 
    a*n[t] - b*p[t],

  tot'[t] == n'[t] + p'[t],

  p[0] == p0,
  n[0] == n0,
  tot[0] == p0 + n0}

sol[gn_, gp_, sn_, sp_, K_, a_, b_, 
   n0_, p0_, tg_, ts_] := 
  First[Evaluate[
    NDSolve[systeqns[gn, gp, sn, sp, 
      K, a, b, n0, p0, tg], {n, p, tot}, {t, 0, tg + ts}, 
     Method -> {"StiffnessSwitching"}, AccuracyGoal -> 1000, 
     WorkingPrecision -> 30, MaxSteps -> Infinity]]];

From which I calculate

avgR[gn_, gp_, sn_, sp_, K_, a_, b_, 
  n0_, p0_, tg_, ts_] := 
 Log[(tot[tg + ts]/tot[0]) /. 
   sol[gn, gp, sn, sp, K, a, b, n0, 
    p0, tg, ts]]

I define some example values

tgex = 5; tsex = 5;gnex = 2;gpex = 2/10;snex = -2;spex = -2/10;bex = 1/10000;n0ex = 1;p0ex = 1/1000;Kex = 10^9;

I am then interested in determining the value of a that maximizes avgR. Graphically, this would be around a=0.2:

vals = Table[{10^pow, 
    avgR[gnex, gpex, snex, spex, Kex, 
     10^pow, bex, n0ex, p0ex, tgex, tsex]}, {pow, -4, 0, 1/10}];
ListLogLinearPlot[vals, PlotRange -> All, Frame -> True, 
 FrameLabel -> {"Switching rate a", "Average population growth rate"},
  PlotStyle -> {Thick}, PlotTheme -> {"Classic", "LargeLabels"}, 
 Joined -> True]

enter image description here

But trying to calculate this value using NMaximize returns an error:

NMaximize[{avgR[gnex, gpex, snex, spex, Kex, 
     a, bex, n0ex, p0ex, tgex, tsex], 
  a > 10^-4 && a < 1}, {a, 1/10}, Method -> "DifferentialEvolution"]
NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0. >>

Any thoughts what I am doing wrong?

Interpolating the values first does work btw, but that is no good for me, since I would eventually like to maximise both a and b, and this would require me calculating too big a grid :

f = Interpolation[vals];
Maximize[{f[a], 10^-4 < a < 1}, a]
{6.507434893938602167476249656, {a -> 
   0.20613498164678882000313750361}}

So I am looking for a solution to get NMaximize working...

The use of a_?NumericQ to restrict evaluation of sol to numeric values of a, which fixes the problem, is explained in the linked resources of this answer to the "pitfalls" question. — Michael E2, Dec 03 '14 at 12:44
@MichaelE2 the evaluation issue seems to be the #1 problem around here. — Yves Klett, Dec 03 '14 at 16:58
@YvesKlett Yes, and it's quite understandable that it is. I wonder if Wolfram would consider mentioning it on documentation pages of the numerics functions. Or perhaps there is a better way to draw users' attention to it. — Michael E2, Dec 03 '14 at 17:05
A bit late, but: usually WhenEvent[] + ParametricNDSolve[] is much more convenient to use in situations like this. — J. M.'s missing motivation, Jul 01 '16 at 01:05

Karsten7 · Accepted Answer · 2014-12-03T07:43:37.640

You have to suppress symbolic evaluations by replacing your definition of avgR with, e.g.

avgR[gn_, gp_, sn_, sp_, K_, a_?NumericQ, b_, n0_, p0_, tg_, ts_] := 
  Log[(tot[tg + ts]/tot[0]) /. sol[gn, gp, sn, sp, K, a, b, n0, p0, tg, ts]]

In NMaximize the 1/10 has to be removed

NMaximize[{avgR[gnex, gpex, snex, spex, Kex, a, bex, n0ex, p0ex, tgex,
tsex], a > 10^-4 && a < 1}, {a}, Method -> "DifferentialEvolution"]

{6.50728, {a -> 0.205164}

If you want to use a starting value and you are only searching for a local maximum you can use

FindMaximum[{avgR[gnex, gpex, snex, spex, Kex, a, bex, n0ex, p0ex, tgex, tsex], 
  a > 10^-4 && a < 1}, {{a, 1/10}}]

{6.50728, {a -> 0.205164}

Using FindMaximum is also much faster than using NMaximize with Method -> "DifferentialEvolution" (6.6 sec vs. 170.1 sec on my PC). NMaximize with Method -> "NelderMead" or Method -> "SimulatedAnnealing" would be faster (13.4 sec and 14.4 sec, respectively) than "DifferentialEvolution". You can get rid of the messages by setting the WorkingPrecision >= the value you set for NDSolve, e.g.

NMaximize[{avgR[gnex, gpex, snex, spex, Kex, a, bex, n0ex, p0ex, tgex, tsex], 
  a > 10^-4 && a < 1}, {a}, Method -> "NelderMead", WorkingPrecision -> 30]

Many thanks for this - really useful that - didn't know about that syntax to suppress symbolic evaluation - great, thx millions! — Tom Wenseleers, Dec 03 '14 at 08:14

Maximising solutions of NDSolve

1 Answers1