如何使用 Python 内置函数 odeint 求解微分方程?

声明:本页面是StackOverFlow热门问题的中英对照翻译,遵循CC BY-SA 4.0协议,如果您需要使用它,必须同样遵循CC BY-SA许可,注明原文地址和作者信息,同时你必须将它归于原作者(不是我):StackOverFlow 原文地址: http://stackoverflow.com/questions/27820725/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me): StackOverFlow

提示:将鼠标放在中文语句上可以显示对应的英文。显示中英文
时间:2020-08-19 02:19:48  来源:igfitidea点击:

How to solve differential equation using Python builtin function odeint?

pythonnumpyscipydifferential-equationsodeint

提问by Physicist

I want to solve this differential equations with the given initial conditions:

我想用给定的初始条件求解这个微分方程:

(3x-1)y''-(3x+2)y'+(6x-8)y=0, y(0)=2, y'(0)=3

the ans should be

y=2*exp(2*x)-x*exp(-x)

答案应该是

y=2*exp(2*x)-x*exp(-x)

here is my code:

这是我的代码:

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()

but what I get is different from the answer. what have I done wrong?

但我得到的与答案不同。我做错了什么?

采纳答案by xnx

There are several things wrong here. Firstly, your equation is apparently

这里有几件事是错误的。首先,你的方程显然是

(3x-1)y''-(3x+2)y'-(6x-8)y=0; y(0)=2, y'(0)=3

(3x-1)y''-(3x+2)y'-(6x-8)y=0;y(0)=2, y'(0)=3

(note the sign of the term in y). For this equation, your analytical solution and definition of y2are correct.

(注意 y 中术语的符号)。对于这个方程,你的解析解和定义y2是正确的。

Secondly, as the @Warren Weckesser says, you must pass 2 parameters as yto g: y[0](y), y[1](y') and return their derivatives, y' and y''.

其次,正如@Warren Weckesser 所说,您必须将 2 个参数传递yg: y[0](y), y[1](y') 并返回它们的导数 y' 和 y''。

Thirdly, your initial conditions are given for x=0, but your x-grid to integrate on starts at -2. From the docs for odeint, this parameter, tin their call signature description:

第三,您的初始条件是为 x=0 给出的,但是您要积分的 x 网格从 -2 开始。从 docs for odeint,这个参数,t在他们的调用签名描述中:

odeint(func, y0, t, args=(),...):

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.

t : array 求解 y 的时间点序列。初始值点应该是这个序列的第一个元素。

So you must integrate starting at 0 or provide initial conditions starting at -2.

因此,您必须从 0 开始积分或提供从 -2 开始的初始条件。

Finally, your range of integration covers a singularity at x=1/3. odeintmay have a bad time here (but apparently doesn't).

最后,您的积分范围涵盖 x=1/3 处的奇点。odeint可能在这里过得不好(但显然没有)。

Here's one approach that seems to work:

这是一种似乎有效的方法:

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()

enter image description here

在此处输入图片说明