How to solve an algebraic equation with InterpolatingFunction?

Question

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?

Answer 1

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]

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. )

Answer 2

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]

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}}*)