python - 求解具有 f (x) 数值但没有解析表达式的颂歌 y'=f (x)

Question

我想用 python 来求解常微分方程，例如y'=f(x)(具有边界条件y(0)=0)。我不知道函数f(x)的解析表达式是什么，相反，我有这个函数的一组点(数据)，用于我想要求解微分方程的域。

我尝试过 odeint 。但当您知道 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) 的值。

Answer 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