Решение оды y'=f (x) с числовыми значениями f (x), но без аналитического выражения
Я хочу численно решить ODE в python, например
y'=f(x) (с граничным условием y(0)=0). Я не знаю, каково аналитическое выражение функции f(x), вместо этого у меня есть набор точек (данных) этой функции для области, в которой я хочу решить дифференциальное уравнение.
Я пробовал с завещание. Но этот метод работает, когда вы знаете явное аналитическое выражение для f(x), что не относится к моему случаю. В частности, у меня есть следующий код с функцией f(x) в массиве (Для простоты я рассматриваю известный f(x), но в моей реальной задаче этот массив f(x) получается из численного моделирования без известного аналитического выражения).
Следующий код не работает, так как odeint думает, что у меня y'=x, а не мои значения f(x):
from scipy.integrate import odeint
import numpy as np
def dy_dx(y, f):
return f #it doesn't work!!
xs = np.linspace(0,10,100)
f = np.sin(xs)*np.exp(-0.1*xs) #data of the function f, but in my real problem I DON'T KNOW THE ANALITICAL SOLUTION! JUST ONLY the points
ys = odeint(dy_dx, 0.0, xs)
В Python должно быть что-то, что может решить эту проблему. По сути, вы решаете оду численно и знаете, каковы значения f(x) в области оды.
И да, функция odeint вызывает dy_dx с аргументами (y,t), а затем, возможно, дополнительными параметрами, указанными в необязательном параметре args=.
Ода, о которой я говорю, является одномерной y(x), xs — это просто массив домена моей оды y'(x)=f(x), поэтому [0,10]. Я думаю, что здесь нечего делать со спаренной системой. Как я уже сказал, я знаю только, каковы значения f(x) в области [0,10], где я хочу решить оду, но я не знаю, что такое аналитическое выражение, поскольку точки f(x ) исходит из другого численного моделирования
Есть ли в другой симуляции вариант плотного вывода? В противном случае вам нужно применить интерполяцию. // Если бы вы рисовали f, как бы вы построили линии?
Да, я могу сделать интерполяцию из вывода f(x). Но все же как в таком случае ввести эту интерполяцию в формализм одеинта?
Ваша задача простая интеграция, почему вы хотите использовать для нее odeint? Кумулятивная квадратурная процедура делает то же самое с меньшими усилиями.
Ответы 1
Вы должны быть в состоянии решить это, используя квадратурные подпрограммы scipy.integrate. Если вы действительно хотите использовать сложную форму, вы должны использовать интерполяцию, например, как в
from scipy.integrate import odeint
from scipy.interpolate import interp1d
import numpy as np
xs = np.linspace(0,10,100+1);
fs = np.sin(xs)*np.exp(-0.1*xs) # = Imag( exp((1j-0.1)*x) )
# the exact anti-derivative of f is
# F = Imag( (exp((1j-0.1)*x)-1)/(1j-0.1) )
# = Imag( -(1j+0.1)*(exp((1j-0.1)*x)-1)/(1.01) )
# = 1/1.01 - exp(-0.1*x)/1.01 * ( cos(x) + 0.1*sin(x) )
def IntF(x): return (1-np.exp(-0.1*x)*(np.cos(x)+0.1*np.sin(x)))/1.01
f = interp1d(xs, fs, kind = "quadratic", fill_value = "extrapolate")
def dy_dx(y, x):
return f(x)
ys = odeint(dy_dx, 0.0, xs)
for x,y in zip(xs, ys): print "%8.4f %20.15f %20.15f"%(x,y,IntF(x))
с первыми 10 строками
x interpolated exact
--------------------------------------------------
0.0000 0.000000000000000 0.000000000000000
0.1000 0.004965420470493 0.004962659238991
0.2000 0.019671988500299 0.019669801188631
0.3000 0.043783730081358 0.043781529336000
0.4000 0.076872788780423 0.076870713937278
0.5000 0.118430993242631 0.118428986914274
0.6000 0.167875357240100 0.167873429717074
0.7000 0.224555718642310 0.224553873611032
0.8000 0.287762489870417 0.287760727322230
0.9000 0.356734939606963 0.356733243391002
1.0000 0.430669760236151 0.430668131955269
Какое измерение имеет
x? В какой форме даны значенияf? Является ли уравнение дляfодновременным дифференциальным уравнением? Если да, пытались ли вы решить ее как связанную систему?