Решение оды y '= f (x) с числовыми значениями f (x), но без аналитического выражения - PullRequest
Новая
       13

Решение оды y '= f (x) с числовыми значениями f (x), но без аналитического выражения

1 голос
/ 08 апреля 2019

Я хочу решить ODE численно в python, например, y '= f (x) (с граничным условием y (0) = 0). Дело в том, что я не знаю, каково аналитическое выражение этой функции f (x), вместо этого у меня есть набор точек (данных) этой функции для области, где я хочу решить дифференциальное уравнение.

Я пробовал с (odeint) https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html. Но этот метод работает, когда вы знаете явное аналитическое выражение для f (x), что не является моим случаем. В частности, у меня есть следующий код с функцией f (x) в массиве (для простоты я считаю, что f (x), НО В МОЕЙ НАСТОЯЩЕЙ ЗАДАЧЕ, ЭТА МОДУЛЬ f (x) СЧИТАЕТСЯ ЧИСЛЕННЫМ МОДЕЛИРОВАНИЕМ, БЕЗ ЗНАКОМ АНАЛИТИЧЕСКИМ EXPRESION).

Я покажу вам код, который я сделал, но не работает, поскольку 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) в области оды. Но я не вижу пути.

Спасибо заранее!

1 Ответ

2 голосов
/ 09 апреля 2019

Вы должны быть в состоянии решить эту проблему, используя квадратурные процедуры 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
...