Нахождение пересечения нуля в питоне - PullRequest
Новая
       34

Нахождение пересечения нуля в питоне

0 голосов
/ 09 мая 2019

Я написал следующий код, чтобы увидеть, в котором t мой ODE "exponential_decay" пересекает нулевую линию.это 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))

. В этом коде мы знаем начальную точку initp = dr[i], мы знаем значение дифференциального уравнения при пересечении нуля desired_answer = odr[i], и мы готовы найти вкоторый t у нас есть этот ответ.Это нормально, и мы получаем ответ по этому коду.ans - это t на пересечении нуля.

Моя проблема в том, что нам делать, если наш ответ ODE, который сейчас равен desired_answer = odr[i], не просто число, а значение в терминах t.

Я имею в виду, используя odr[i], мы читаем файл данных и затем читаем числа.Теперь рассмотрим, что у нас есть что-то вроде odr = 0.1 * t, 0.12 *t, 0.23 *t etc., которое не является числом и является функцией с точки зрения t.

1 Ответ

1 голос
/ 10 мая 2019

Это не самое эффективное использование интерфейса solve_ivp.Вы можете получить свой результат автоматически, используя механизм event. 

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]

Даже если вы хотите использовать свой собственный решатель (или вызов решателя), вы можете получить точную кусочно-полиномиальную интерполяцию для конкретного метода, используяопция плотного вывода, 

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

Чтобы получить корни функции, зависящей от времени, вам просто нужно закодировать эту зависимость от времени, например, как в 

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

или 

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

Если вы хотите получить достоверные результаты, вы должны явно контролировать допуски atol, rtol.

...