如何在 scipy.integrate.ode （或其他）函数中编写初始条件？

Question

我正在尝试使用 python scipy.integrate.odeint 函数求解微分方程并将其与 mathcad 解进行比较。

所以我的方程是 u'' + 0.106u'+ 0.006u = 0, 我遇到的问题是初始条件 u(0)=0 and u'(1)=1。我不明白如何设置 u'(1)=1。

Python代码：

from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np

def eq(u,t):
    return [u[1], -0.106*u[1]-0.006*u[0]] #return u' and u''

time = np.linspace(0, 10) 
u0 = [0,1] # initial conditions
Z = odeint(eq,u0,time) </code>


plt.plot(time, Z)
plt.xticks(range(0,10))
plt.grid(True)
plt.xlabel('time')
plt.ylabel('u(t)')
plt.show()

Mathcad 代码：

u''(t) + 0.106*u'(t) +0.006*u(t) = 0
u(0) = 0
u'(1) = 1
u := Odesolve(t,10)

Mathcad 图表如下所示：

https://pp.userapi.com/c850032/v850032634/108079/He1JsQonhpk.jpg
这是标准具。

我的 python 输出是：

https://pp.userapi.com/c850032/v850032634/10809c/KB_HDekc8Fk.jpg
看起来确实相似，但显然 u(t) 不正确。

Answer 1

这是一道边值题，需要用到solve_bvp

from scipy.integrate import solve_bvp, solve_ivp
import matplotlib.pyplot as plt
import numpy as np

def eq(t,u): return [u[1], -0.106*u[1]-0.006*u[0]] #return u' and u''
def bc(u0,u1): return [u0[0], u1[1]-1 ]

res = solve_bvp(eq, bc, [0,1], [[0,1],[1,1]], tol=1e-3)
print res.message

# plot the piecewise polynomial interpolation, 
# using the last segment for extrapolation
time = np.linspace(0, 10, 501) 
plt.plot(time, res.sol(time)[0], '-', lw=2, label="BVP extrapolated")
# solve the extended solution as IVP solution on [1,10]
ivp = solve_ivp(eq, time[[0,-1]], res.y[:,0], t_eval=time)
plt.plot(time, ivp.y[0], '-', lw=2, label="IVP from BVP IC")

# plot decorations
plt.xticks(range(0,11))
plt.grid(True)
plt.legend()
plt.xlabel('time')
plt.ylabel('u(t)')
plt.show()

请注意，延续是通过从给定区间 [0,1] 到 [0,10] 的外推法得出的，并且 1 处的值具有 1e-3 的公差。因此，通过使用 solve_ivp 并将 t=1 处的计算值作为初始值，可以在大区间内获得更好的结果。此示例中的差异约为 0.01.

如何在 scipy.integrate.ode （或其他）函数中编写初始条件？

How to write initial conditions in scipy.integrate.ode (or another) function?

python

scipy

nonlinear-optimization

differential-equations

odeint

Python代码：

Mathcad 代码：

Mathcad 图表如下所示：

我的 python 输出是：