1
$\begingroup$

I'm trying to fit three differential equations to two data sets, using the code developed by @Cesareo, as follows:

tk = {0, 3, 6, 9, 15, 18, 21, 24};
y1 = {0.974, 0.492, 0.829, 0.554, 0.208, 0.138, 0.199, 0.0893};
y2 = {0.915094, 0.736097, 0.793694, 0.833664, 1, 0.99578, 0.897964, 0.214499};


fitness[gamma_?NumberQ, c_?NumberQ, beta_?NumberQ, b_?NumberQ, r_?NumberQ, n_?NumberQ, m_?NumberQ, p0_?NumberQ, T0_?NumberQ] := 
Module[{p, f, T, odes, ODES, solode, cinits, lambda = 1, t}, 
odes = {p'[t] == 0.47 p[t] (1 - p[t]/5.29) - gamma p[t] f[t], f'[t] == -c f[t] + beta p[t] (1 - f[t]), T'[t] == b (1 - T[t]) (1 - 1/(1 + (r/f[t])^n)) - m (1/(1 + (r/f[t])^n))}; 
cinits = {p[0] == p0, f[0] == 0.001, T[0] == T0}; 
ODES = Join[odes, cinits]; 
solode = Quiet@NDSolve[ODES, {p, f, T}, {t, tk[[1]], tk[[8]]}][[1]];
Return[Sum[Norm[{y1[[k]], lambda y2[[k]]} - Evaluate[{p[t], lambda T[t]} /. solode /. {t -> tk[[k]]}]], {k, 1, 8}]]]


restrs = {0 < gamma < 1 && 0 < c < 1 && 0 < beta < 1 && 0 < b < 1 && 0 < r < 1 && 1 < n < 15 && 0 < m < 1 && 0.8 < p0 < 1 && 0.8 < T0 < 1.5};
vars = {gamma, c, beta, b, r, n, m, p0, T0};


sol = Quiet@FindMinimum[Join[{fitness[gamma, c, beta, b, r, n, m, p0, T0]}, restrs], {{gamma, 0.01}, {c, 0.01}, {beta, 0.01}, {b, 0.01}, {r, 0.01}, {n, 2}, {m, 0.01}, {p0, 0.85}, {T0, 0.85}}, 
StepMonitor :> Print[{gamma, c, beta, b, r, n, m, p0, T0, fitness[gamma, c, beta, b, r, n, m, p0, T0]}], MaxIterations -> 10, WorkingPrecision -> 30]


odes = {p'[t] == 0.47 p[t] (1 - p[t]/5.29) - gamma p[t] f[t], f'[t] == -c f[t] + beta p[t] (1 - f[t]), T'[t] == b (1 - T[t]) (1 - 1/(1 + (r/f[t])^n)) - m (1/(1 + (r/f[t])^n))} /. sol[[2]];
cinits = {p[0] == p0, f[0] == 0.001, T[0] == T0} /. sol[[2]];
ODES = Join[odes, cinits];
solode = NDSolve[ODES, {p, f, T}, {t, tk[[1]], tk[[8]]}][[1]];


pl1 = Plot[Evaluate[p[t] /. solode], {t, tk[[1]], tk[[8]]}, PlotStyle -> Blue, AxesOrigin -> {0, 0}];
pl2 = ListPlot[Transpose[{tk, y1}], PlotStyle -> Blue, AxesOrigin -> {0, 0}];
Show[pl1, pl2]


pl3 = Plot[Evaluate[f[t] /. solode], {t, tk[[1]], tk[[8]]}, PlotStyle -> Green, AxesOrigin -> {0, 0}]


pl4 = Plot[Evaluate[T[t] /. solode], {t, tk[[1]], tk[[8]]}, PlotStyle -> Red, AxesOrigin -> {0, 0}];
pl5 = ListPlot[Transpose[{tk, y2}], PlotStyle -> Red, AxesOrigin -> {0, 0}];
Show[pl4, pl5]

My problem is that, I really don't know the restrictions on parameters as they are specified in restrs in the code. How can one include the built-in Manipulate for the parameters, in order to vary them and observe at the same time their effects on the fits?

$\endgroup$

1 Answer 1

1
$\begingroup$

You could add slider for the lower and upper limit of each entry in restrs. Since you have 9 of these, then you need 18 sliders

Mathematica graphics

Note code does not check that lower limit is less than upper limit. You should add such a check in your code at one point or make sure the lower limit slider stays below the upper limit slides. I also added ContinuousAction -> False, since I found the code slow. This means Manipulate will update when you release your hand from the slider and not as the slider is moving. You are free to remove this if you want.

Code

tk = {0, 3, 6, 9, 15, 18, 21, 24}; y1 = {0.974, 0.492, 0.829, 0.554, 
  0.208, 0.138, 0.199, 0.0893}; y2 = {0.915094, 0.736097, 0.793694, 
  0.833664, 1, 0.99578, 0.897964, 0.214499};
fitness[gamma_?NumberQ, c_?NumberQ, beta_?NumberQ, b_?NumberQ, 
   r_?NumberQ, n_?NumberQ, m_?NumberQ, p0_?NumberQ, T0_?NumberQ] := 
  Module[{p, f, T, odes, ODES, solode, cinits, lambda = 1, t},
   odes = {p'[t] == 0.47  p[t]  (1 - p[t]/5.29) - gamma  p[t]  f[t], 
     f'[t] == -c  f[t] + beta  p[t]  (1 - f[t]), 
     T'[t] == 
      b  (1 - T[t])  (1 - 1/(1 + (r/f[t])^n)) - 
       m  (1/(1 + (r/f[t])^n))};
   cinits = {p[0] == p0, f[0] == 0.001, T[0] == T0};
   ODES = Join[odes, cinits];
   solode = Quiet@NDSolve[ODES, {p, f, T}, {t, tk[[1]], tk[[8]]}][[1]];
   Return[
    Sum[Norm[{y1[[k]], lambda  y2[[k]]} - 
       Evaluate[{p[t], lambda  T[t]} /. 
          solode /. {t -> tk[[k]]}]], {k, 1, 8}]]];

Add Manipulate

Manipulate[
 Module[{restrs, gamma, c, beta, b, r, n, m, p0, T0, vars, sol, odes, 
   p, T, t, f, cinits, ODES, solode, pl1, pl2, pl3, pl4, pl5, p1, 
   p2},
  
  vars = {gamma, c, beta, b, r, n, m, p0, T0};
  restrs = {gammaLower < gamma < gammaUpper && cLower < c < cUpper && 
     betaLower < beta < betaUpper && bLower < b < bUpper && 
     rLower < r < rUpper && nLower < n < nUpper && 
     mLower < m < mUpper && p0Lower < p0 < p0Upper && 
     T0Lower < T0 < T0Upper};
  sol = Quiet@
    FindMinimum[
     Join[{fitness[gamma, c, beta, b, r, n, m, p0, T0]}, 
      restrs], {{gamma, 0.01}, {c, 0.01}, {beta, 0.01}, {b, 0.01}, {r,
        0.01}, {n, 2}, {m, 0.01}, {p0, 0.85}, {T0, 0.85}}, 
     MaxIterations -> 10, WorkingPrecision -> 30];
  odes = {p'[t] == 0.47  p[t]  (1 - p[t]/5.29) - gamma  p[t]  f[t], 
     f'[t] == -c  f[t] + beta  p[t]  (1 - f[t]), 
     T'[t] == 
      b  (1 - T[t])  (1 - 1/(1 + (r/f[t])^n)) - 
       m  (1/(1 + (r/f[t])^n))} /. sol[[2]];
  cinits = {p[0] == p0, f[0] == 0.001, T[0] == T0} /. sol[[2]];
  ODES = Join[odes, cinits];
  solode = NDSolve[ODES, {p, f, T}, {t, tk[[1]], tk[[8]]}][[1]];
  pl1 = Plot[Evaluate[p[t] /. solode], {t, tk[[1]], tk[[8]]}, 
    PlotStyle -> Blue, AxesOrigin -> {0, 0}, ImageSize -> 400];
  pl2 = ListPlot[Transpose[{tk, y1}], PlotStyle -> Blue, 
    AxesOrigin -> {0, 0}, ImageSize -> 400];
  p1 = Show[pl1, pl2];
  pl3 = Plot[Evaluate[f[t] /. solode], {t, tk[[1]], tk[[8]]}, 
    PlotStyle -> Green, AxesOrigin -> {0, 0}, ImageSize -> 400]; 
  pl4 = Plot[Evaluate[T[t] /. solode], {t, tk[[1]], tk[[8]]}, 
    PlotStyle -> Red, AxesOrigin -> {0, 0}, ImageSize -> 400]; 
  pl5 = ListPlot[Transpose[{tk, y2}], PlotStyle -> Red, 
    AxesOrigin -> {0, 0}, ImageSize -> 400]; p2 = Show[pl4, pl5];
  Row[{p1, p2, pl3}]
  ],
 {{gammaLower, 0, "lower gamma"}, -10, 10, 1, Appearance -> "Labeled"},
 {{gammaUpper, 1, "upper gamma"}, -10, 10, 1, Appearance -> "Labeled"},
 {{cLower, 0, "lower c"}, -10, 10, 1, Appearance -> "Labeled"},
 {{cUpper, 1, "upper c"}, -10, 10, 1, Appearance -> "Labeled"},
 {{betaLower, 0, "lower beta"}, -10, 10, 1, Appearance -> "Labeled"},
 {{betaUpper, 1, "upper beta"}, -10, 10, 1, Appearance -> "Labeled"},
 {{bLower, 0, "lower b"}, -10, 10, 1, Appearance -> "Labeled"},
 {{bUpper, 1, "upper b"}, -10, 10, 1, Appearance -> "Labeled"},
 {{rLower, 0, "lower r"}, -10, 10, 1, Appearance -> "Labeled"},
 {{rUpper, 1, "upper r"}, -10, 10, 1, Appearance -> "Labeled"},
 {{nLower, 1, "lower n"}, -10, 20, 1, Appearance -> "Labeled"},
 {{nUpper, 15, "upper n"}, -10, 20, 1, Appearance -> "Labeled"},
 {{mLower, 0, "lower m"}, -10, 20, 1, Appearance -> "Labeled"},
 {{mUpper, 1, "upper m"}, -10, 20, 1, Appearance -> "Labeled"},
 {{p0Lower, 0.8, "lower p0"}, -5, 5, .1, Appearance -> "Labeled"},
 {{p0Upper, 1, "upper p0"}, -5, 5, .1, Appearance -> "Labeled"},
 {{T0Lower, 0.8, "lower T0"}, -5, 5, .1, Appearance -> "Labeled"},
 {{T0Upper, 1.5, "upper T0"}, -5, 5, .1, Appearance -> "Labeled"},
 ContinuousAction -> False,
 TrackedSymbols :> {gammaLower, gammaUpper, cLower, cUpper, betaLower,
    betaUpper, bLower, bLower, rLower, rUpper, nLower, nUpper, mLower,
    mUpper, p0Lower, p0Upper, T0Lower, T0Upper}
 ]
$\endgroup$
1
  • 1
    $\begingroup$ @Anovice You just need to move the slides to set the lower limit and the upper limit for each of your 9 variables. But as I said, there is no check now in code to detect that lower limit is smaller than upper limit automatically. This is something that can be added later. $\endgroup$
    – Nasser
    Commented Feb 2 at 10:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.