TransWikia.com

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

Mathematica Asked on July 19, 2021

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.

boundary conditionsdifferential equationsnumerical integration

2 Answers

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, 
    "StartingInitialConditions" -> Thread[Through[dvars@0] == params]},
  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" -> 
    Thread[Through[dvars@0] == params /. icsNM]},
 WorkingPrecision -> MachinePrecision]

Answered by Michael E2 on July 19, 2021

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]}

Figure 1 Unfortunately there is no any method to compare with. Also maximal error is of $3times 10^{-3}$ and mean error is of $8times10^{-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"}] ]

Figure 2

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}]]

Figure 3

Answered by Alex Trounev on July 19, 2021

Add your own answers!

Ask a Question

Get help from others!

Recent Answers

Recent Questions

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP