这里对使用python求解常微分方程提供两种思路,一种是自己编程实现欧拉法,改进欧拉法或者四阶龙格库塔,这样有助于理解上述三种数值计算方法的原理;一种是调用python已有的库,不再重复造轮子。
本文对上述两种思路都给出代码示例,并进行比较;同时针对单个微分方程和含有多个微分方程的微分方程组给出代码示例
1. 常微分方程定义
凡含有参数,未知函数和未知函数导数 (或微分) 的方程,称为微分方程。
未知函数是一元函数的微分方程称作常微分方程
未知数是多元函数的微分方程称作偏微分方程。
微分方程中出现的未知函数最高阶导数的阶数,称为微分方程的阶数。
2. 调用现有的库
scipy中提供了用于解常微分方程的函数odeint(),完整的调用形式如下:
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)
实际使用中,还是主要使用前三个参数,即微分方程的描写函数、初值和需要求解函数值对应的的时间点。接收数组形式。这个函数,要求微分方程必须化为标准形式,即
,实际操作中会发现,高阶方程的标准化工作,其实是解微分方程最主要的工作。
示例1:单个微分方程
import math
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def func(y, t):
return t * math.sqrt(y)
YS=odeint(func,y0=1,t=np.arange(0