Python 中求解微分方程的问题

Question

我正在努力完成this，这里我要解决五个普通的问题使用 odeint 的微分方程并重现该任务中给出的数字。

这是我的代码：

import scipy as sp
import scipy.interpolate as ip
import numpy as np
import matplotlib.pyplot as pl

d = 8.64
Mu1 = 4.95*10**2
Mu2 = 4.95*10**(-2)
vs = 0.12
vd = 1.23
w = 10**(-3)
k1 = 2.19*10**(-4)
k2 = 6.12*10**(-5)
k3 = 0.997148
k4 = 6.79*10**(-2)

p0 = 1.00
sigmas0 = 2.01
sigmad0 = 2.23
alphas0 = 2.20
alphad0 = 2.26

hs = (sigmas0-(sigmas0**(2)-k3*alphas0*(2*sigmas0-alphas0))**(1/2))/k3
cs = (alphas0-hs)/2
ps = k4*(hs**2)/cs

t_points = [ 1000, 1850, 1950, 1980, 2000, 2050, 2080, 2100, 2120, 2150, 2225, 2300, 2500, 5000 ]
y_points = [ 0.0, 0.0, 1.0, 4.0, 5.0, 8.0, 10.0, 10.5, 10.0, 8.0, 3.5, 2.0, 0.0, 0.0 ]

t1 = np.array(t_points)
y1 = np.array(y_points)

new_length = 1000
new_t = np.linspace(t1.min(), t1.max(), new_length)
new_y2 = ip.pchip_interpolate(t1, y1, new_t)

pl.plot(t_points,y_points,'o', new_t,new_y2)
pl.show()

ft = sp.interpolate.interp1d(new_t, new_y2)

def equations(x, t1):

    p = x[0]
    alphad = x[1]
    alphas = x[2]
    sigmad = x[3] 
    sigmas = x[4]

    dpdt = (ps-p)/d + ft/Mu1
    dalphaddt = (1/vd)*(k2-w*(alphad-alphas))
    dalphasdt = (1/vs)*(w*(alphad-alphas)-k2)
    dsigmaddt = (1/vd)*(k1-w*(sigmad-sigmas))
    dsigmasdt = (1/vs)*(w*(sigmad-sigmas)-k1-(ps-p)/d*Mu2)

    return [dpdt, dalphaddt, dalphasdt, dsigmaddt, dsigmasdt]

solve =  sp.integrate.odeint(equations, [p0, alphad0, alphas0, sigmad0, sigmas0], t1)

好像是这部分：

dpdt = (ps-p)/d + ft/Mi1

错误，我不知道如何解决。

错误说：

TypeError: unsupported operand type(s) for /: 'interp1d' and 'float'.

非常感谢任何想法和建议。

编辑：当我应用 meowgoesthedog 建议的更改时，出现错误：

Traceback (most recent call last):

  File "<ipython-input-5-324757833872>", line 1, in <module>
    runfile('E:/Data/Project 2/project2.py', wdir='E:/Data/Project 2')

  File "D:\Programs\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 668, in runfile
    execfile(filename, namespace)

  File "D:\Programs\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 108, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "E:/Data/Project 2/project2.py", line 59, in <module>
    solve =  odeint(equations, [p0, alphad0, alphas0, sigmad0, sigmas0], t1)

  File "D:\Programs\Anaconda3\lib\site-packages\scipy\integrate\odepack.py", line 233, in odeint
    int(bool(tfirst)))

  File "E:/Data/Project 2/project2.py", line 51, in equations
    dpdt = (ps-p)/d + ft(t1)/Mu1

  File "D:\Programs\Anaconda3\lib\site-packages\scipy\interpolate\polyint.py", line 79, in __call__
    y = self._evaluate(x)

  File "D:\Programs\Anaconda3\lib\site-packages\scipy\interpolate\interpolate.py", line 664, in _evaluate
    below_bounds, above_bounds = self._check_bounds(x_new)

  File "D:\Programs\Anaconda3\lib\site-packages\scipy\interpolate\interpolate.py", line 696, in _check_bounds
    raise ValueError("A value in x_new is above the interpolation "

ValueError: A value in x_new is above the interpolation range.

`

Answer 1

根据interp1d's documentation：

ynew = f(xnew) # use interpolation function returned by interp1d

它returns 一个函数/可调用对象，其取值x 和returns f(x) 的内插值。在你的情况下 "x" = t:

dpdt = (ps-p)/d + ft(t1)/Mu1   # pass t1 to ft to obtain interpolated value

更新

此新错误是由于 odeint 在 t 超出的值处对函数 f(t) 进行采样t_points 的值。这是纠正错误所必需的，并且没有选项可以阻止 odeint 这样做。但是，我们可以 extrapolate f(t) 超出提供的样本，使用 InterpolatedUnivariateSpline:
```
from scipy.interpolate import InterpolatedUnivariateSpline

...    

ft = InterpolatedUnivariateSpline(t1, y1, k=1)
```
与 interp1d 一样，此 returns 具有相同签名的函数。但是，应用此修复后，结果变为：

这当然是不正确的。

您已在函数外部声明 hs, cs, ps 为常量。事实上，它们是 alpha* 和 sigma* 变量的函数，因此必须在每次调用 equation:

def equations(x, t):

    p = x[0]
    alphad = x[1]
    alphas = x[2]
    sigmad = x[3]
    sigmas = x[4]

    hs = (sigmas-(sigmas**(2)-k3*alphas*(2*sigmas-alphas))**(1/2))/k3
    cs = (alphas-hs)/2
    ps = k4*(hs**2)/cs

    dpdt = (ps-p)/d + ft(t)/Mu1
    dalphaddt = (1/vd)*(k2-w*(alphad-alphas))
    dalphasdt = (1/vs)*(w*(alphad-alphas)-k2)
    dsigmaddt = (1/vd)*(k1-w*(sigmad-sigmas))
    dsigmasdt = (1/vs)*(w*(sigmad-sigmas)-k1-(ps-p)/d*Mu2)

    return [dpdt, dalphaddt, dalphasdt, dsigmaddt, dsigmasdt]

结果现在与练习中的图表相匹配...几乎。

您将 t1 作为水平轴变量传递给了 odeint。它只有 14 个元素，对于平滑输出来说太少了。改为传递 new_t：
```
solve = ig.odeint(equations, [p0, alphad0, alphas0, sigmad0, sigmas0], new_t)
```
现在的结果完全符合预期！

Python 中求解微分方程的问题

Troubles with solving differental equations in Python

python

class

scipy

odeint