# Is Shooting the best NDSolve method for two point conditions, and how to improve its accuracy?

Hi I have a two point conditions problem (controlled SIR epidemics via Pontryagin BVP), which is supposed to depend heavily on initial conditions, which break the problem into several cases (choosing one case is the uncomprehensible part of my program in the beginning). I realized that even my best NDSolve Shooting method may get wrong answers (for example for the initial values, see printout below); in fact, from the comments below it may be that using NDSolve in two-point problems may not be ideal.

In light of the difficulties mentioned also in the answers, I am going to add two more aspects to the quest. The optimization problem also requires to flatten the green i curve" (infections) by a state constraint i \leq iM, which may be handled by adding WhenEvent[i[t] > iM, i[t] -> iM-.0001] to the equations.

    Clear["Global*"]; \[Alpha] = 5/10; \[Rho] = 137/100; \[Gamma] =
1/18; R0 = 28/10; \[Beta] = \[Gamma] R0;
bl = 4 \[Beta]/(10); Rl =
bl/\[Gamma]; h = ((# - 1 - Log[#])/#) &; g = (1 - 1/#) &;
g0 = g[R0]; h0 = h[R0];
gs = g[Rl]; hs = h[Rl]; iM2 = (4 gs + hs)/5;(*CHOOSE iM*)iM = gs;
s0f = (Exp[# (1 - iM) - 1]/#) &;
s0 = N[s0f[Rl]] + .07; i0 = 1 - s0 // N;
tf = 80;
pw = Re[(#/\[Rho])^(1/(-1 + \[Rho]))] &;
arg = Min[\[Beta],
Max[\[Beta] - pw[#], bl]] &(*doesn't work with piecewise*);
upon[t_] := arg[i[t] s[t] (pi[t] - ps[t])];
sys = {odes,
bcs} = {{s'[t] == -u[t]*s[t]*i[t],
i'[t] == i[t] (u[t]*s[t] - \[Gamma]),
ps'[t] == u[t] i[t] (ps[t] - pi[t]),
pi'[t] == u[t] i[t] (ps[t] - pi[t]) - \[Alpha] + pi[t] \[Gamma],
WhenEvent[i[t] > iM, i[t] -> iM-.0001]}, {i[0] == i0, s[0] == s0,
ps[tf] == 0, pi[tf] == \[Alpha]/\[Gamma]}};
dvars = {i, s, ps, pi};
sol = NDSolve[sys /. u[t] -> upon[t], dvars, {t, 0, tf},
Method -> "Shooting", WorkingPrecision -> MachinePrecision];
Print["Rl=",
Rl // N, " hard i_M inter= ", {h[Rl], g[Rl]} //
N, " s0=", s0, " >I=",
1/R0 // N, " i0=", (i /. First[sol])[0], "=", i0, " iM=", iM // N]
uop = Plot[{Evaluate[{(upon[t]/\[Beta]), bl/\[Beta], i[t], iM,
s[t](*,\[Psi][t]*)} /. sol]}, {t, 0, tf},
PlotLegends -> {"control =upon[t]/\[Beta]",
"lower bound for control", "i[t]=infectious",
"iM=max allowable for infection", "s[t]=susceptibles"(*,
"\[Psi][t] (determines control"*)},
Epilog -> {Black, Dashed, Line[{{0, 1}, {tf, 1}}]},
PlotLabel ->
Style[Framed["Plot of control upon/\[Beta]"], 16, Blue,
Background -> Lighter[Yellow]], AxesLabel -> Automatic]


The optimality is supposed to come from the choice upon[t_] := arg[i[t] s[t] (pi[t] - ps[t])] see https://elischolar.library.yale.edu/cgi/viewcontent.cgi?article=1213&context=cowles-discussion-paper-series

The goal of the lowerbound bl' is to pinpoint the optimal period for total confinement, when the control is a lower bang". A sharp transition to the lower bang is sure to happen when $$\rho$$ becomes 1, but seems difficult to get numerically.

Even with bigger $$\rho$$, bad things happen: note that the computed i[0] in my example is twice the one specified by the initial condition.

The plots, assuming tf is big enough, are supposed to show an increasing then decreasing i(infection). Sometimes https://arxiv.org/abs/2007.00318 there may be infectious with two optimal peaks, and possibly several peaks for non-optimal arbitrary government policies might be possible.

• Briefly: (1) There might be no solution. (2) If there is a solution, you need "StartingInitialConditions" that are good enough; problems in which there is little knowledge of what good starting points could be called "difficult." (3) Method -> "FiniteElement" is now an option for nonlinear BVPs (in V12.?); it can be used either to solve the system or get good starting points. (4) ParametricNDSolve[] can be used for greater control over shooting (one could use NMinimize on the norm of the residual of the boundary conditions instead of FindRoot). (5) There's the chasing method. Feb 13, 2021 at 20:22
• @florin NDSolve is not right method for optimization problem. Please, have a look my answer to this post mathematica.stackexchange.com/questions/239572/… Feb 17, 2021 at 15:19

Here's a way that would work for ODE systems that do not develop singularities, assuming the desired BVP has a solution. We use NMinimize to hunt down decent starting initial conditions.

sys = {odes, bcs} = {{s'[t] == -u[t]*s[t]*i[t],
i'[t] == i[t] (u[t]*s[t] - γ),
ps'[t] == u[t] i[t] (ps[t] - pi[t]),
pi'[t] == u[t] i[t] (ps[t] - pi[t]) - α + pi[t] γ},
{i[0] == i0, s[0] == s0, ps[tf] == 0, pi[tf] == (α/γ)}};
dvars = {i, s, ps, pi};

params = {i1, s1, ps1, pi1};
psol = ParametricNDSolveValue[
{odes, Through[dvars@0] == params} /. u[t] -> upon[t],
dvars, {t, 0, tf}, params,
Method -> {"Shooting", Method -> Automatic,
WorkingPrecision -> MachinePrecision]

bcerr[i1_?NumericQ, s1_?NumericQ, ps1_?NumericQ, pi1_?NumericQ] :=
Total[(
bcs /. Equal -> Subtract /.
Thread[dvars -> psol @@ {i1, s1, ps1, pi1}]
)^2];
{norm, icsNM} =
NMinimize[{bcerr @@ params , And @@ Thread[params > 0]},
params,
WorkingPrecision -> MachinePrecision, PrecisionGoal -> 6,
MaxIterations -> 10]
(*  {2.28422, {i1 -> 1.13849, s1 -> 0.000107863, ps1 -> 12.6763, pi1 -> 13.2513}}  *)

sol = NDSolve[
sys /. u[t] -> upon[t],
dvars, {t, 0, tf},
Method -> {"Shooting", Method -> "ExplicitRungeKutta",
"StartingInitialConditions" ->
WorkingPrecision -> MachinePrecision]

• How we know that result of this code is numerical optimal solution of the SIR epidemic model with time-dependent control? Probably we need some reference to compare. Feb 17, 2021 at 15:09
• @AlexTrounev I suppose you're referring to BVPs to possibly having multiple solutions and choosing one that optimizes some measure of utility. I don't know anything about "optimal" SIR solutions. What I did doesn't try to optimize anything, just find one solution to the BVP. (NMinimize is used to cast a broader net to try to locate a solution, it being less dependent on having a good starting position than FindRoot.) The final solution is the result of the built-in "Shooting" method, which gives no error in this case and so should be as reliable as NDSolve. Feb 17, 2021 at 15:52
• Any way (+1) for the original solution. Just pay attention, it is not working for $\rho = 137/100$ and around. And I am not able to check how NDSolve[] accepts u[t]. There are no messages, but may be NDSolve[] just ignores this function. Feb 17, 2021 at 20:58
• @AlexTrounev Thanks. It seems to be a "difficult" BVP. NMinimize struggles to minimize the norm of the residuals of the BVP in the given problem, and never could meet the precision goal IIRC. I limited the max iterations to 10 because more didn't help much and used a lot of time. I gave up on an automatic method that would solve for all parameters. I'm not surprised that adjustments would be necessary for different ranges of parameter values. Perhaps it is just too sensitive. I don't feel I have time for an in-depth analysis of it. Feb 17, 2021 at 22:04
• @AlexTrounev thanks for pointing out NDSolve may not be ideal method for BVP, and for your advice, which I'm studying now. For basics of UNCONTROLLED SIR, see en.wikipedia.org/wiki/… The controlled problem, giving rise to Pontryagin type problems is a recent research topic from last year -- see for example arxiv.org/abs/2007.00318 and it has been tackled sometimes via the package bocop.org/download I am trying to do it with Mathematica Feb 18, 2021 at 17:26

This is kind of optimization problem, and it can be solved with NMinimize[] using colocation method and Haar wavelets as follows. Define parameters (we change parameter $$\rho$$ to some reliable value)

tf = 80; \[Alpha] = 5/10; \[Rho] =
137/100; \[Gamma] = 1/18; R0 =
28/10; \[Beta] = \[Gamma] \
R0;
bl = 4 \[Beta]/(10); Rl =
bl/\[Gamma]; H = ((# - 1 - Log[#])/#) &; g = (1 - 1/#) &;
g0 = g[R0]; h0 = H[R0];
gs = g[Rl]; hs = H[Rl]; iM2 = (gs + hs)/2;
iM = iM2;
s0f = (Exp[# (1 - iM) - 1]/#) &;
s0 = N[s0f[Rl]] + .05;
Print["Rl=",
Rl // N, " hard i_M inter= ", {H[Rl], g[Rl]} //
N, " s0=", s0, " >I=", 1/R0 // N, " i0=", i0 = 1 - s0 // N, " iM=",
iM // N]



Define collocation points, wavelets, functions and derivatives

A = 0; B = 1; J = 4; M = 2^J; dx = (B - A)/(2*M);
h1[x_] := Piecewise[{{1, A <= x <= B}, {0, True}}];
p1[x_, n_] := (1/n!)*(x - A)^n;
h[x_, k_, m_] := Piecewise[{{1, Inequality[k/m, LessEqual, x, Less,
(1 + 2*k)/(2*m)]}, {-1, Inequality[(1 + 2*k)/(2*m), LessEqual, x, Less,
(1 + k)/m]}}, 0]
p[x_, k_, m_, n_] := Piecewise[{{0, x < k/m}, {(-(k/m) + x)^n/n!,
Inequality[k/m, LessEqual, x, Less, (1 + 2*k)/(2*m)]},
{((-(k/m) + x)^n - 2*(-((1 + 2*k)/(2*m)) + x)^n)/n!,
(1 + 2*k)/(2*m) <= x <= (1 + k)/m},
{((-(k/m) + x)^n + (-((1 + k)/m) + x)^n - 2*(-((1 + 2*k)/(2*m)) + x)^n)/n!,
x > (1 + k)/m}}, 0]
xl = Table[A + l*dx, {l, 0, 2*M}]; xcol = Table[(xl[[l - 1]] + xl[[l]])/2,
{l, 2, 2*M + 1}];
iH1[x_] := Sum[aiH[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
a0*h1[x];
iH0[x_] := Sum[aiH[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
a0*p1[x, 1] + a1;
sH1[x_] := Sum[asH[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
b0*h1[x];
sH0[x_] := Sum[asH[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
b0*p1[x, 1] + b1;
psH1[x_] := Sum[apsH[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
c0*h1[x];
psH0[x_] := Sum[apsH[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
c0*p1[x, 1] + c1;
piH1[x_] := Sum[apiH[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
d0*h1[x];
piH0[x_] := Sum[apiH[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
d0*p1[x, 1] + d1;


Define variables, function u[t], constrains and equations

var = Flatten[
Table[{aiH[i, j], asH[i, j], apiH[i, j], apsH[i, j]}, {j, 0, J,
1}, {i, 0, 2^j - 1, 1}]]; varX =
Join[{a0, a1, b0, b1, c0, c1, d0, d1}, var];
pw[x_] := Re[(x/\[Rho])^(1/(-1 + \[Rho]))];
arg[x_] := Min[\[Beta], Max[\[Beta] - pw[x], bl]];
u[t_] := arg[iH0[t] sH0[t] (piH0[t] - psH0[t])];
cons0 = {iH0[0] - i0 == 0, sH0[0] - s0 == 0, psH0[1] == 0,
piH0[1] - (\[Alpha]/\[Gamma]) == 0};
cons1 = Flatten[
Table[{iH0[t] >= 0, sH0[t] >= 0, piH0[t] >= 0, psH0[t] >= 0}, {t,
xcol}]];
cons = Join[cons0, cons1];

eqn = {-s'[t]/tf - u[t]*s[t]*i[t], -i'[t]/tf +
i[t] (u[t]*s[t] - \[Gamma]), -ps'[t]/tf +
u[t] i[t] (ps[t] - pi[t]), -pi'[t]/tf +
u[t] i[t] (ps[t] - pi[t]) - \[Alpha] + pi[t] \[Gamma]} /. {i ->
iH0, i' -> iH1, s -> sH0, s' -> sH1, ps -> psH0, ps' -> psH1,
pi -> piH0, pi' -> piH1};
eq =  Flatten[Table[eqn, {t, xcol}]];


Solve minimization problem

solM = NMinimize[{Norm[eq], cons}, varX]

(*{0.0104057, {a0 -> 0.0954815, a1 -> 0.00508882, b0 -> -0.375913,
b1 -> 0.994911, c0 -> -3.90017, c1 -> 3.90017, d0 -> -0.677472,
d1 -> 9.67747, aiH[0, 0] -> 0.0660761, asH[0, 0] -> -0.063544,
apiH[0, 0] -> 1.00814, apsH[0, 0] -> 0.405986,
aiH[0, 1] -> 0.188154, asH[0, 1] -> 0.0376435,
apiH[0, 1] -> 0.487586, apsH[0, 1] -> 0.593654,
aiH[1, 1] -> 0.0444829, asH[1, 1] -> -0.0189795,
apiH[1, 1] -> 1.0024, apsH[1, 1] -> 0.150677, aiH[0, 2] -> 0.06544,
asH[0, 2] -> 0.11632, apiH[0, 2] -> 0.725894, apsH[0, 2] -> 0.87757,
aiH[1, 2] -> -0.0234203, asH[1, 2] -> -0.00566019,
apiH[1, 2] -> 0.338734, apsH[1, 2] -> 0.30369,
aiH[2, 2] -> 0.0252932, asH[2, 2] -> -0.0104613,
apiH[2, 2] -> 0.258855, apsH[2, 2] -> 0.0964121,
aiH[3, 2] -> -0.0265804, asH[3, 2] -> -0.00918152,
apiH[3, 2] -> 0.756419, apsH[3, 2] -> 0.0304156,
aiH[0, 3] -> -0.0745347, asH[0, 3] -> 0.137549,
apiH[0, 3] -> 0.650927, apsH[0, 3] -> 0.875131,
aiH[1, 3] -> 0.170202, asH[1, 3] -> -0.0411005,
apiH[1, 3] -> -0.257863, apsH[1, 3] -> -0.248828,
aiH[2, 3] -> -0.0160979, asH[2, 3] -> 0.00409634,
apiH[2, 3] -> 0.190301, apsH[2, 3] -> 0.219103,
aiH[3, 3] -> -0.0198636, asH[3, 3] -> -0.0159907,
apiH[3, 3] -> 0.172247, apsH[3, 3] -> 0.103688,
aiH[4, 3] -> 0.0401597, asH[4, 3] -> -0.0287653,
apiH[4, 3] -> -0.169186, apsH[4, 3] -> -0.225655,
aiH[5, 3] -> 0.0146636, asH[5, 3] -> 0.00353028,
apiH[5, 3] -> 0.260637, apsH[5, 3] -> 0.124237,
aiH[6, 3] -> 0.000467526, asH[6, 3] -> -0.00932397,
apiH[6, 3] -> 0.24602, apsH[6, 3] -> -0.0247289,
aiH[7, 3] -> 0.0109868, asH[7, 3] -> 0.00531618,
apiH[7, 3] -> 0.533452, apsH[7, 3] -> 0.0562497,
aiH[0, 4] -> -0.0161501, asH[0, 4] -> 0.0650699,
apiH[0, 4] -> 0.2376, apsH[0, 4] -> 0.391809,
aiH[1, 4] -> -0.0423136, asH[1, 4] -> 0.0705367,
apiH[1, 4] -> 0.412222, apsH[1, 4] -> 0.478024,
aiH[2, 4] -> 0.0484061, asH[2, 4] -> -0.0108476,
apiH[2, 4] -> -0.00221506, apsH[2, 4] -> -0.00930998,
aiH[3, 4] -> 0.0791673, asH[3, 4] -> -0.0115428,
apiH[3, 4] -> -0.11338, apsH[3, 4] -> -0.0820615,
aiH[4, 4] -> 0.0275533, asH[4, 4] -> 0.0115671,
apiH[4, 4] -> 0.102044, apsH[4, 4] -> 0.128047,
aiH[5, 4] -> -0.0124024, asH[5, 4] -> -0.00230487,
apiH[5, 4] -> 0.0969294, apsH[5, 4] -> 0.0958279,
aiH[6, 4] -> -0.00120586, asH[6, 4] -> -0.00779982,
apiH[6, 4] -> 0.0704392, apsH[6, 4] -> 0.0433793,
aiH[7, 4] -> -0.0307616, asH[7, 4] -> -0.0103363,
apiH[7, 4] -> 0.10174, apsH[7, 4] -> 0.0517215,
aiH[8, 4] -> -0.000430565, asH[8, 4] -> -0.0314585,
apiH[8, 4] -> -0.153533, apsH[8, 4] -> -0.199178,
aiH[9, 4] -> 0.00501252, asH[9, 4] -> 0.00718202,
apiH[9, 4] -> 0.0937051, apsH[9, 4] -> 0.0750744,
aiH[10, 4] -> -0.00331805, asH[10, 4] -> 0.00282372,
apiH[10, 4] -> 0.132389, apsH[10, 4] -> 0.0815513,
aiH[11, 4] -> 0.020478, asH[11, 4] -> -0.000615826,
apiH[11, 4] -> 0.132089, apsH[11, 4] -> 0.0428073,
aiH[12, 4] -> 0.0103351, asH[12, 4] -> -0.00509096,
apiH[12, 4] -> 0.104502, apsH[12, 4] -> -0.0156213,
aiH[13, 4] -> -0.0103834, asH[13, 4] -> -0.0039744,
apiH[13, 4] -> 0.135814, apsH[13, 4] -> -0.0174982,
aiH[14, 4] -> -0.0082928, asH[14, 4] -> 0.00002658,
apiH[14, 4] -> 0.219744, apsH[14, 4] -> 0.0166508,
aiH[15, 4] -> 0.0374625, asH[15, 4] -> 0.0114181,
apiH[15, 4] -> 0.311253, apsH[15, 4] -> 0.0345916}}*)



Plot numerical solution and check errors for every colocation point

{Plot[
Evaluate[{iH0[t/tf], sH0[t/tf], psH0[t/tf], piH0[t/tf]} /.
solM[[2]]], {t, 0, tf},
PlotLegends -> {"i, NMinimize", "s, NMinimize", "ps, NMinimize",
"pi, NMinimize"}, AxesOrigin -> {0, 0}],

Plot[Evaluate[{iH0[t/tf], sH0[t/tf]} /. solF], {t, 0, tf},
PlotLegends -> {"i, FindRoot", "s, FindRoot"}, AxesOrigin -> {0, 0},
PlotStyle -> {Blue, Red}, PlotPoints -> 200]}


Unfortunately there is no any method to compare with. Also maximal error is of $$3\times 10^{-3}$$ and mean error is of $$8\times10^{-5}$$, and it corresponds to $$(\Delta x)^2=1/32^2$$, nevertheless we need some reference data. For this we run FindRoot[] with initial point at solM as

solF = FindRoot[Join[Table[eq[[i]] == 0, {i, Length[eq]}], cons0],
Table[{varX[[i]], varX[[i]] /. solM[[2]]}, {i, Length[varX]}],
MaxIterations -> 1000];


Finally we plot 2 solutions and check improvement

Show[Plot[Evaluate[{iH0[t/tf], sH0[t/tf]} /. solF], {t, 0, tf},
PlotLegends -> {"i, FindRoot", "s, FindRoot"}, AxesOrigin -> {0, 0},
PlotStyle -> {Blue, Red}, PlotPoints -> 200], Plot[Evaluate[{iH0[t/tf], sH0[t/tf]} /. solM[[2]]], {t, 0, tf},
PlotLegends -> {"i, NMinimize", "s, NMinimize"}] ]


Show[Plot[Evaluate[{piH0[t/tf], psH0[t/tf]} /. solF], {t, 0, tf},
PlotLegends -> {"pi, FindRoot", "ps, FindRoot"},
AxesOrigin -> {0, 0}, PlotStyle -> {Blue, Red}, PlotPoints -> 200],
Plot[Evaluate[{piH0[t/tf], psH0[t/tf]} /. solM[[2]]], {t, 0, tf},
PlotLegends -> {"pi, NMinimize", "ps, NMinimize"},
AxesOrigin -> {0, 0}]]


• (+1) Perhaps you meant to format the Plot code at the end? Feb 18, 2021 at 0:21
• @MichaelE2 Thank you! Next step is to use FindRoot[] with initial point at solM. Feb 18, 2021 at 10:48