用 f (x) 的数值求解一个 ode y'=f (x) 但没有分析表达式 [英] Solving an ode y'=f (x) with numerical values of f (x) but without analitical expresion

问题描述

我想在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将 numpy 导入为 npdef dy_dx(y, f):return f #它不起作用！！xs = np.linspace(0,10,100)f = np.sin(xs)*np.exp(-0.1*xs) #函数 f 的数据，但在我真正的问题中我不知道解析解！只有积分ys = odeint(dy_dx, 0.0, xs)

Python 中一定有什么东西可以解决这个问题.基本上，您正在以数字方式求解 ode，并且您知道 ode 域中 f(x) 的值.

解决方案

您应该能够使用 scipy.integrate 的正交例程来解决这个问题.如果你真的想使用复杂的形式，你必须使用插值，例如

from scipy.integrate import odeint从 scipy.interpolate 导入 interp1d将 numpy 导入为 npxs = np.linspace(0,10,100+1);fs = np.sin(xs)*np.exp(-0.1*xs) # = Imag( exp((1j-0.1)*x) )# f 的精确反导数是# F = Imag( (exp((1j-0.1)*x)-1)/(1j-0.1) )# = 图像( -(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.01f = interp1d(xs, fs, kind="quadratic", fill_value="extrapolate")def dy_dx(y, x):返回 f(x)ys = odeint(dy_dx, 0.0, xs)对于 zip(xs, ys) 中的 x,y:打印 "%8.4f %20.15f %20.15f"%(x,y,IntF(x))

前 10 行

 x 内插精确--------------------------------------------------0.0000 0.0000000000000000 0.0000000000000000.1000 0.004965420470493 0.0049626592389910.2000 0.019671988500299 0.0196698011886310.3000 0.043783730081358 0.0437815293360000.4000 0.076872788780423 0.0768707139372780.5000 0.118430993242631 0.1184289869142740.6000 0.167875357240100 0.1678734297170740.7000 0.224555718642310 0.2245538736110320.8000 0.287762489870417 0.2877607273222300.9000 0.356734939606963 0.3567332433910021.0000 0.430669760236151 0.430668131955269

I want to solve an ODE numerically in python, like y'=f(x) (with boundary condition y(0)=0). I don't know what the analytical expression of function f(x), instead I have a set of points (data) of this function for the domain where I want to solve the differential equation.

I have tried with odeint. But this method works when you know the explicit analytical expression for f(x), which is not my case. In particular, I have the following code with the function f(x) in an array (For simplicity I consider a known f(x), but in my real problem this array f(x) comes from a numerical simulation without known analytical expression).

The following code doesn't work, since odeint thinks that I have y'=x, and not my values 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)

There must be something in Python that can solve this. Basically you are solving the ode numerically and you know what the values of f(x) in the domain of the ode.

解决方案

You should be able to solve this using the quadrature routines of scipy.integrate. If you really want to use the complicated form, you have to use interpolation, for instance as in

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))

with the first 10 lines

    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

这篇关于用 f (x) 的数值求解一个 ode y'=f (x) 但没有分析表达式的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！