如何使用 Python 内置函数 odeint 求解微分方程?
How to solve differential equation using Python builtin function odeint?
我想用给定的初始条件求解这个微分方程:
(3x-1)y''-(3x+2)y'+(6x-8)y=0, y(0)=2, y'(0)=3
答案应该是
y=2*exp(2*x)-x*exp(-x)
这是我的代码:
def g(y,x):
y0 = y[0]
y1 = y[1]
y2 = (6*x-8)*y0/(3*x-1)+(3*x+2)*y1/(3*x-1)
return [y1,y2]
init = [2.0, 3.0]
x=np.linspace(-2,2,100)
sol=spi.odeint(g,init,x)
plt.plot(x,sol[:,0])
plt.show()
但是我得到的和答案不一样。
我做错了什么?
这里有几处错误。首先,你的等式显然是
(3x-1)y''-(3x+2)y'-(6x-8)y=0; y(0)=2, y'(0)=3
(注意 y 中项的符号)。对于这个方程,你的解析解和 y2
的定义是正确的。
其次,正如@Warren Weckesser 所说,您必须将 2 个参数作为 y
传递给 g
:y[0]
(y)、y[1]
(y')和 return 它们的导数,y' 和 y''。
第三,你的初始条件是 x=0,但你的 x-grid 积分从 -2 开始。来自 odeint
的文档,这个参数,t
在他们的调用签名描述中:
odeint(func, y0, t, args=(),...)
:
t : array
A sequence of time points for which to solve for y. The initial
value point should be the first element of this sequence.
因此您必须从 0 开始积分或提供从 -2 开始的初始条件。
最后,您的积分范围涵盖了 x=1/3 处的奇点。 odeint
在这里可能过得不好(但显然没有)。
这是一种似乎有效的方法:
import numpy as np
import scipy as sp
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def g(y, x):
y0 = y[0]
y1 = y[1]
y2 = ((3*x+2)*y1 + (6*x-8)*y0)/(3*x-1)
return y1, y2
# Initial conditions on y, y' at x=0
init = 2.0, 3.0
# First integrate from 0 to 2
x = np.linspace(0,2,100)
sol=odeint(g, init, x)
# Then integrate from 0 to -2
plt.plot(x, sol[:,0], color='b')
x = np.linspace(0,-2,100)
sol=odeint(g, init, x)
plt.plot(x, sol[:,0], color='b')
# The analytical answer in red dots
exact_x = np.linspace(-2,2,10)
exact_y = 2*np.exp(2*exact_x)-exact_x*np.exp(-exact_x)
plt.plot(exact_x,exact_y, 'o', color='r', label='exact')
plt.legend()
plt.show()
我想用给定的初始条件求解这个微分方程:
(3x-1)y''-(3x+2)y'+(6x-8)y=0, y(0)=2, y'(0)=3
答案应该是
y=2*exp(2*x)-x*exp(-x)
这是我的代码:
def g(y,x):
y0 = y[0]
y1 = y[1]
y2 = (6*x-8)*y0/(3*x-1)+(3*x+2)*y1/(3*x-1)
return [y1,y2]
init = [2.0, 3.0]
x=np.linspace(-2,2,100)
sol=spi.odeint(g,init,x)
plt.plot(x,sol[:,0])
plt.show()
但是我得到的和答案不一样。 我做错了什么?
这里有几处错误。首先,你的等式显然是
(3x-1)y''-(3x+2)y'-(6x-8)y=0; y(0)=2, y'(0)=3
(注意 y 中项的符号)。对于这个方程,你的解析解和 y2
的定义是正确的。
其次,正如@Warren Weckesser 所说,您必须将 2 个参数作为 y
传递给 g
:y[0]
(y)、y[1]
(y')和 return 它们的导数,y' 和 y''。
第三,你的初始条件是 x=0,但你的 x-grid 积分从 -2 开始。来自 odeint
的文档,这个参数,t
在他们的调用签名描述中:
odeint(func, y0, t, args=(),...)
:
t : array A sequence of time points for which to solve for y. The initial value point should be the first element of this sequence.
因此您必须从 0 开始积分或提供从 -2 开始的初始条件。
最后,您的积分范围涵盖了 x=1/3 处的奇点。 odeint
在这里可能过得不好(但显然没有)。
这是一种似乎有效的方法:
import numpy as np
import scipy as sp
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def g(y, x):
y0 = y[0]
y1 = y[1]
y2 = ((3*x+2)*y1 + (6*x-8)*y0)/(3*x-1)
return y1, y2
# Initial conditions on y, y' at x=0
init = 2.0, 3.0
# First integrate from 0 to 2
x = np.linspace(0,2,100)
sol=odeint(g, init, x)
# Then integrate from 0 to -2
plt.plot(x, sol[:,0], color='b')
x = np.linspace(0,-2,100)
sol=odeint(g, init, x)
plt.plot(x, sol[:,0], color='b')
# The analytical answer in red dots
exact_x = np.linspace(-2,2,10)
exact_y = 2*np.exp(2*exact_x)-exact_x*np.exp(-exact_x)
plt.plot(exact_x,exact_y, 'o', color='r', label='exact')
plt.legend()
plt.show()