求解带有 f (x) 数值但没有解析表达式的 ode y'=f (x)答案

【问题标题】：Solving an ode y'=f (x) with numerical values of f (x) but without analitical expresion求解带有 f (x) 数值但没有解析表达式的 ode y'=f (x)
【发布时间】：2019-08-29 16:28:27
【问题描述】：

我想在 python 中以数值方式求解 ODE，例如 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 中一定有一些东西可以解决这个问题。基本上你是在用数字求解 ode 并且你知道 ode 域中 f(x) 的值是什么。

【问题讨论】：

x 有什么维度？ f 的值以什么形式给出？ f 的方程是联立微分方程吗？如果是，您是否尝试将其作为耦合系统来解决？
是的，odeint 函数调用dy_dx 并带有参数(y,t)，然后可能在args= 可选参数中给出了附加参数。
我说的ode是一维y(x)，xs就是我的ode y'(x)=f(x)的域的数组，所以[0,10] .我认为这里与耦合系统无关。正如我所说，我只知道要解决 ode 的域 [0,10] 中 f(x) 的值是多少，但我不知道分析表达式是什么，因为点 f(x ) 来自其他数值模拟
其他模拟有密集输出选项吗？否则，您需要应用插值。 // 如果你要绘制f，你将如何生成线条？
是的，我可以从输出 f(x) 进行插值。但是，在这种情况下，如何在 odeint 的形式中引入这种插值？

标签： python ode

【解决方案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

【讨论】：