用 scipy 在 python 中求解二维微分方程答案

【问题标题】：solving two dimension-differential equations in python with scipy用 scipy 在 python 中求解二维微分方程
【发布时间】：2016-04-09 13:52:40
【问题描述】：

我是 python 的新手。我有一个简单的微分系统，它由两个变量和两个微分方程和初始条件x0=1, y0=2：

dx/dt=6*y
dy/dt=(2t-3x)/4y

现在我正在尝试解决这两个微分方程，我选择odeint。这是我的代码：

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

def func(z,b):
    x, y=z
    return [6*y, (b-3*x)/(4*y)]    

z0=[1,2]
t = np.linspace(0,10,11)
b=2*t
xx=odeint(func, z0, b)
pl.figure(1)
pl.plot(t, xx[:,0])
pl.legend()
pl.show()

但是结果不正确，有错误提示：

Excess work done on this call (perhaps wrong Dfun type).
Run with full_output = 1 to get quantitative information.

我不知道我的代码有什么问题以及如何解决它。任何帮助都会对我有用。

【问题讨论】：

与stackoverflow.com/questions/34616775/… 中的错误消息和上下文相同。只是一个更简单的功能。
stackoverflow.com/questions/28092456/… 中的可能答案。这谈到了相同的lsoda 警告，我正在运行您的功能。他说问题不在于代码，而在于 ode 的行为。
谢谢你，我正在阅读链接的问题。我想知道，我是否使用 odeint 权利来解决这两个微分方程？我从未使用 odeint 来求解两个微分方程
你使用了正确的接口，但是你没有实现公布的ODE。解决y'(t)=f(2t,y(t)) 和y'(t)=f(t,y(t)) 存在差异，这不是通过将样本点加倍来解决的。
integrate.odeint(func, [1.,2.],np.linspace(0,.6,10)) 运行良好。但是当我将范围提高到 0.7 时，我开始收到此警告。正如LutzL 所说，y 太接近于 0，而你的导数正在爆炸。

标签： python scipy ode

【解决方案1】：

应用技巧对除以 y 的除法进行去奇异化，打印所有 ODE 函数评估，绘制两个分量，并使用修改后的代码使用正确的微分方程

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

def func(z,t):
    x, y=z
    print t,z
    return [6*y, (2*t-3*x)*y/(4*y**2+1e-12)]    

z0=[1,2]
t = np.linspace(0,1,501)
xx=odeint(func, z0, t)
pl.figure(1)
pl.plot(t, xx[:,0],t,xx[:,1])
pl.legend()
pl.show()

您会看到在t=0.64230232515 处假设了y=0 的奇点，其中y 在其顶点处表现得像平方根函数。没有办法越过那个奇点，因为y 的斜率趋于无穷大。此时解不再连续可微，因此这是解的极值点。常量延续是去奇异化的产物，不是有效的解决方案。

【讨论】：