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

Question

我是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.

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

Answer 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， where 的y行为类似于其顶点处的平方根函数。没有办法越过那个奇点，因为 的斜率y趋于无穷大。此时，解不再连续可微，因此这是解的极值点。持续延续是去单一化的产物，而不是有效的解决方案。