1
$\begingroup$

How to solve the equation f[t]==0 for the function f[t] defined by differential equation,i.e. if f[t] is no analytic? That means that f[t] is the form InterpolatingFunction

For a dynamic system defined by differential equations, I find the inflection points, extremes et.c. For instance, let the system is the SIR model.

 β = 1.5; γ = 1/3; i0 = 0.001;
  sol = NDSolve[
  {s'[t] == -β*s[t]*i[t],
  i'[t] == β*s[t]*i[t] - γ*i[t],
  r'[t] == γ*i[t],
  i[0] == i0, r[0] == 0, s[0] == 1 - i0},
 {s, i, r},
 {t, 0, 104}]

The system has no exact algebraic solution. The output is in the form

{{s -> InterpolatingFunction[...], i -> InterpolatingFunction[...], r -> InterpolatingFunction[...]}

However. the Solve don't work as so as NSolve, Reduce et.c. Some failed attempts to find numerically a maxima of i[t], i.e. solving the equation $$ \beta i(t) s(t)-\gamma i(t)=0. $$

NSolve[((β*s[t]*i[t] - γ*i[t]) /. sol) == 0., t]
NSolve[(Evaluate@(β*s[t]*i[t] - γ*i[t]) /. sol ) == 0., t]
NSolve[(Evaluate@(β*s[t]*i[t] - γ*i[t]) /. sol /. 
t -> τ ) == 0., τ]
... 

How to solve it?

$\endgroup$
3
$\begingroup$

Yes, NSolve (and Solve) can't handle this equation involving InterpolatingFunction (at least now), but it's not hard to circumvent. The basic approach is to observe the plot to find the rough position of the root and use FindRoot:

eq = (β s[t] i[t] - γ i[t])

Plot[eq /. sol // Evaluate, {t, 0, 104}, PlotRange -> All]

enter image description here

FindRoot[eq /. sol, {t, 0, 15}]

(* {t -> 7.13011} *)

Notice the initial guess of FindRoot needs to be set carefully. You may also want to try the functions in this post.

Then here comes the slightly advanced solution:

Last@Reap@NDSolve[{s'[t] == -β s[t] i[t], 
    i'[t] == β s[t] i[t] - γ i[t], r'[t] == γ i[t], i[0] == i0, 
    r[0] == 0, s[0] == 1 - i0, WhenEvent[eq == 0 // Evaluate, Sow@t]}, {s, i, r}, {t, 0, 
    104}, WorkingPrecision -> 32]

(* {{7.1301355405901063248884664090806}} *)

The WorkingPrecision option is necessary or fake root will be found. (As we can see from the plot, the solution is so close to zero when t gets large. )

| improve this answer | |
$\endgroup$
  • 2
    $\begingroup$ If only maxima of i[t] are required, WhenEvent[i'[t] < 0, Sow@t] works nicely (no extra WorkingPrecision needed). $\endgroup$ – Chris K Apr 10 at 12:57
  • $\begingroup$ To the downvoter, I am interested in what's wrong with my answer, would you please elaborate. I'm not trying to complain here, I just want to improve my answer if possible. $\endgroup$ – xzczd Apr 11 at 4:01
  • $\begingroup$ Event location is indeed the best method here. $\endgroup$ – J. M.'s technical difficulties May 19 at 7:42
1
$\begingroup$

It's a difficult task!

Look at the plot

Plot[((\[Beta]*s[t]*i[t] - \[Gamma]*i[t]) /. sol)[[1]], {t, 0, 104},PlotRange -> All]

enter image description here

There "zeros" of your equation aren't well defined.

Perhaps you better look for the minimal time your equation approaches zero?

NMinimize[{t, ((\[Beta]*s[t]*i[t] - \[Gamma]*i[t]) /.  sol)[[1]]^2 < .001,10< t < 104}, t]
(*{13.2505, {t -> 13.2505}}*)
| improve this answer | |
$\endgroup$
  • $\begingroup$ However, the NSolve[(Evaluate@(\[Beta]*s[t]*i[t] - \[Gamma]*i[t]) /. sol ) == 0., t, 3.] does not work as well. I think the problem is not in the actual function `i[t]' but its presentation as 'InterpolatingFunction' $\endgroup$ – Slepecky Mamut Apr 10 at 11:54
  • 1
    $\begingroup$ Perhaps interpolation function are the problem for NSolve but you can solve it with FindRoot or NMinimize. Look at the plot. Perhaps you have to help Mathematica a little bit. Are you looking for zeros t>20 or 5<t<10? $\endgroup$ – Ulrich Neumann Apr 10 at 12:04
  • 1
    $\begingroup$ The InterpolatingFunctions are the problem. A simple example: f = FunctionInterpolation[x, {x, -6, 6}]; g = FunctionInterpolation[2 x, {x, -6, 6}]; NSolve[f[x] + g[x] == 0, x] Please see my answer for more details. $\endgroup$ – xzczd Apr 10 at 12:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.