# How to solve an algebraic equation with InterpolatingFunction?

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 == i0, r == 0, s == 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?

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 == i0,
r == 0, s == 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. )

• If only maxima of i[t] are required, WhenEvent[i'[t] < 0, Sow@t] works nicely (no extra WorkingPrecision needed). – Chris K Apr 10 '20 at 12:57
• 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. – xzczd Apr 11 '20 at 4:01
• Event location is indeed the best method here. – J. M.'s torpor May 19 '20 at 7:42

Look at the plot

Plot[((\[Beta]*s[t]*i[t] - \[Gamma]*i[t]) /. sol)[], {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)[]^2 < .001,10< t < 104}, t]
(*{13.2505, {t -> 13.2505}}*)

• 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' – Slepecky Mamut Apr 10 '20 at 11:54
• 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? – Ulrich Neumann Apr 10 '20 at 12:04
• 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. – xzczd Apr 10 '20 at 12:14