Python_learning4：python中的解常微分方程 odeint函数 - [转载自Python科学计算]

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

版权声明：此文非博主原创，以超链接形式标明文章原始出处和作者信息及本声明：

链接：点击打开转载源

Link：http://www.blogbus.com/blankdesktop-logs/87973246.html

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

常用形式
odeint(func, y0, t,args,Dfun)
一般这种形式就够用了。
下面是官方的例子，求解的是
D(D(y1))-t*y1=0
为了方便，采取D=d/dt。如果我们令初值
y1(0) = 1.0/3**(2.0/3.0)/gamma(2.0/3.0)
D(y1)(0) = -1.0/3**(1.0/3.0)/gamma(1.0/3.0)
这个微分方程的解y1=airy(t)。
 
令D(y1)=y0，就有这个常微分方程组。
D(y0)=t*y1
D(y1)=y0
 
python求解该微分方程。
>>> from scipy.integrate import odeint
>>> from scipy.special import gamma, airy
>>> y1_0 = 1.0/3**(2.0/3.0)/gamma(2.0/3.0)
>>> y0_0 = -1.0/3**(1.0/3.0)/gamma(1.0/3.0)
>>> y0 = [y0_0, y1_0]
>>> def func(y, t):
...     return [t*y[1],y[0]]
>>> def gradient(y,t):
...     return [[0,t],[1,0]]
>>> x = arange(0,4.0, 0.01)
>>> t = x
>>> ychk = airy(x)[0]
>>> y = odeint(func, y0, t)
>>> y2 = odeint(func, y0, t, Dfun=gradient)
>>> print ychk[:36:6]
[ 0.355028  0.339511  0.324068  0.308763  0.293658  0.278806]
>>> print y[:36:6,1]
[ 0.355028  0.339511  0.324067  0.308763  0.293658  0.278806]
>>> print y2[:36:6,1]
[ 0.355028  0.339511  0.324067  0.308763  0.293658  0.278806]
 
得到的解与精确值相比，误差相当小。
=======================================================================================================
 
args是额外的参数。
用法请参看下面的例子。这是一个洛仑兹曲线的求解，并且用matplotlib绘出空间曲线图。（来自《python科学计算》）
from scipy.integrate import odeint
import numpy as np
def lorenz(w, t, p, r, b):
    # 给出位置矢量w，和三个参数p, r, b 计算出
    # dx/dt, dy/dt, dz/dt 的值
    x, y, z = w
    # 直接与lorenz 的计算公式对应
    return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])
t = np.arange(0, 30, 0.01) # 创建时间点
# 调用ode 对lorenz 进行求解, 用两个不同的初始值
track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))
track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0))
# 绘图
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(track1[:,0], track1[:,1], track1[:,2])
ax.plot(track2[:,0], track2[:,1], track2[:,2])
plt.show()
=========================================================================== 
scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)
计算常微分方程（组）
使用 FORTRAN库odepack中的lsoda解常微分方程。这个函数一般求解初值问题。
 
参数：

• func : callable(y, t0, ...)    计算y在t0 处的导数。
• y0 : 数组    y的初值条件（可以是矢量）
• t : 数组    为求出y，这是一个时间点的序列。初值点应该是这个序列的第一个元素。
• args : 元组    func的额外参数
• Dfun : callable(y, t0, ...)    函数的梯度(Jacobian)。即雅可比多项式。
• col_deriv : boolean.    True，Dfun定义列向导数（更快），否则Dfun会定义横排导数
• full_output : boolean    可选输出，如果为True 则返回一个字典，作为第二输出。
• printmessg : boolean    是否打印convergence 消息。

 
返回： y : array, shape (len(y0), len(t))
数组，包含y值，每一个对应于时间序列中的t。初值y0 在第一排。
infodict : 字典，只有full_output == True 时，才会返回。
字典包含额为的输出信息。
键值：

• ‘hu’   vector of step sizes successfully used for each time step.
• ‘tcur’ vector with the value of t reached for each time step. (will always be at least as large as the input times).
• ‘tolsf’ vector of tolerance scale factors, greater than 1.0, computed when a request for too much accuracy was detected.
• ‘tsw’  value of t at the time of the last method switch (given for each time step)
• ‘nst’   cumulative number of time steps
• ‘nfe’   cumulative number of function evaluations for each time step
• ‘nje’   cumulative number of jacobian evaluations for each time step
• ‘nqu’  a vector of method orders for each successful step.
• ‘imxer’index of the component of largest magnitude in the weighted local error vector (e / ewt) on an error return, -1 otherwise.
• ‘lenrw’  the length of the double work array required.
• ‘leniw’  the length of integer work array required.
• ‘mused’a vector of method indicators for each successful time step: 1: adams (nonstiff), 2: bdf (stiff)

其他参数，官方网站和文档都没有明确说明。相关的资料，暂时也找不到。