1

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?

xzczd
  • 65,995
  • 9
  • 163
  • 468
Slepecky Mamut
  • 1,762
  • 7
  • 15

2 Answers2

3

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

xzczd
  • 65,995
  • 9
  • 163
  • 468
1

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}}*)
Ulrich Neumann
  • 53,729
  • 2
  • 23
  • 55
  • 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
  • 1
    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
  • 1
    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