Solving Interpolating Function

Questions
, ,
1

Hi,

I’ve solved a system of three differential equations to produce three interpolating functions. I need to find the value of t when the first interpolating function is equal to 1. Any idea’s on how to do this? Thanks.

ScreenShot2014-04-17at11.55.17.png
ScreenShot2014-04-17at11.55.17.png1234×820 76.7 KB
Tha

2

If you copy your code into a “spikey” windows instead of posting a screen shot, other people can copy it.

3 
A = 0.01
B = 0.00001
eq1 = h'[t] == -A *h[t]* z[t]*r[t]
eq2 = z'[t] == A *h[t]*z[t] - B *h[t]*z[t]
eq3 = r'[t] == B *h[t] *z[t]
% sol = First[
  h /. NDSolve[{eq1, eq2, eq3, h[0] == 1000000, z[0] == 1, 
     r[0] == 0}, {h, r, z}, {t, 0, 20}]]
sol = First[
  h /. NDSolve[{eq1, eq2, eq3, h[0] == 99.9, z[0] == 0.1, 
     r[0] == 0}, {h, r, z}, {t, 0, 20}]]
sol2 = First[
  z /. NDSolve[{eq1, eq2, eq3, h[0] == 99.9, z[0] == 0.1, 
     r[0] == 0}, {h, r, z}, {t, 0, 20}]]

I put the text above in and got three interpolating functions out. I need to solve for t when the first interpolating function is equal to 0.

4 
A = 0.01;
B = 0.00001;
eq1 = h'[t] == -A*h[t]*z[t]*r[t];
eq2 = z'[t] == A*h[t]*z[t] - B*h[t]*z[t];
eq3 = r'[t] == B*h[t]*z[t];
sol = First[
   h /. NDSolve[{eq1, eq2, eq3, h[0] == 1000000, z[0] == 1, 
      r[0] == 0}, {h, r, z}, {t, 0, 20}]];


FindRoot[sol[t] == 0, {t, 10^-3}]
{t -> 0.00288353}

There seem to be a number of points where the solution is zero, looking at a graph of the log of the solution.

5

The multiple zero crossings are just an artifact of using machine precision. It looks like h is monotonically decreasing and 0 is a horizontal asymptote that it can get arbitrarily close to, so one just needs to decide at what tolerance h will be considered to vanish.

In[1]:= A = 1/100;
B = 1/100000;
eq1 = h'[t] == -A*h[t]*z[t]*r[t];
eq2 = z'[t] == A*h[t]*z[t] - B*h[t]*z[t];
eq3 = r'[t] == B*h[t]*z[t];
sol = NDSolveValue[{eq1, eq2, eq3, h[0] == 1000000, z[0] == 1, 
    r[0] == 0}, h, {t, 0, 20}, WorkingPrecision -> 200, 
   Method -> "StiffnessSwitching"];

In[7]:= sol /@ {3/1000, 4/1000, 5/1000} // N

Out[7]= {1.85407*10^-11, 3.89607*10^-20, 8.19073*10^-29}
6

Oh, now I see the original question was about crossing 1, not 0.

In[10]:= FindRoot[sol[t] == 1, {t, 0, 20}, WorkingPrecision -> 200] // N

Out[10]= {t -> 0.00176324}

In[11]:= sol[t /. %]

Out[11]= 1.

Alternatively, could use WhenEvent:

In[12]:= Reap[NDSolveValue[{eq1, eq2, eq3, h[0] == 1000000, z[0] == 1, 
      r[0] == 0, WhenEvent[h[t] == 1, Sow[t]]}, h, {t, 0, 20}, 
     WorkingPrecision -> 200, Method -> "StiffnessSwitching"]][[2, 1]] // N

Out[12]= {0.00176324}