0

I have written the following code to see in which t my ODE "exponential_decay" crosses the zero line. this is Brent Method.

odr, hr, dr, cr, m = np.genfromtxt('data.txt',unpack=True)
n=0
with open('RDE_nob_trans.txt', 'w') as d:
  for i in range(len(dr)):
      c = cr[i]
      initp = dr[i]
      exponential_decay = lambda t, y: -(1/(1+t)) * (2 *(1-y)* (-2 + (y/c)) + 3 - 3 * y)
      t_span = [0, 1] # Interval of integration
      y0 = [initp] # Initial state: y(t=t_span[0])=2
      desired_answer = odr[i]
      sol_ode = solve_ivp(exponential_decay, t_span, y0) # IVP solution

      f_sol_ode = interp1d(sol_ode.t, sol_ode.y) # Build interpolated function
      ans = brentq(lambda x: f_sol_ode(x) - desired_answer, t_span[0], t_span[1])
      d.write("{0} {1} {2} {3} {4}\n".format(hr[i], dr[i], cr[i], m[i], ans))

In this code we know the initial point initp = dr[i], we know the value of the differential equation at the zero crossing desired_answer = odr[i], and we are willing to find in which t we have this answer. It is OK and we get the answer by this code. ans is the t at the zero crossing.

My problem is what should we do if our answer of ODE which now is desired_answer = odr[i] is not just a number and is a value in terms of t.

I mean using odr[i] we are reading the data file and subsequently reading numbers. Now consider we have someting like odr = 0.1 * t, 0.12 *t, 0.23 *t etc. which is not a number and is a function in terms of t.

Ma Y
  • 696
  • 8
  • 19

1 Answers1

1

This is not the most efficient use of the solve_ivp interface. You can get your result automatically using the event mechanism.

def event(t,y): return y-desired_answer;
event.terminal = True;

sol_ode = solve_ivp(exponential_decay, t_span, y0, events=event) # IVP solution
ans = sol_ode.t[-1]

Even if you want to use your own solver (or solver call), you can get a method-specific and accurate piecewise polynomial interpolation using the dense output option,

sol_ode = solve_ivp(exponential_decay, t_span, y0, dense_output=True) # IVP solution
f_sol_ode = sol_ode.sol

To get the roots of a function that depends on the time, you just have to code that time dependence, for instance as in

def event(t,y): return y-0.23*t;

or

ans = brentq(lambda t: f_sol_ode(t) - 0.23*t, t_span[0], t_span[1])

If you want reliable results, you should explicitly control the tolerances atol, rtol.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • the upper part of your code does not work for me, full of error. I appreciate if you write a complete code, because python is taking about no effect on `event` solver etc. about the lower part. it works, but a problem that I have , I think i forgot to mention in the question is: the code just work for `t_span = [0, 1]` and I cannot extend the range of that for example `-1,1` or `0,2` etc. in some values it says `ValueError: f(a) and f(b) must have different signs` – Ma Y May 10 '19 at 08:34
  • Can you give two lines of example parameters? I corrected the syntax errors in the parameter name and boolean constant. And one can probably not take the first component of the sol function, only of its value arrays. – Lutz Lehmann May 10 '19 at 08:54
  • As to the error, you can apply a bracketing method like bisection, regula falsi, Dekker, Müller or Brent only to intervals where you have observed a sign change in the target function. The event mechanism does this automatically, checking for sign changes during an internal step, you have to do this explicitly. – Lutz Lehmann May 10 '19 at 08:58