自变量只有一个的微分方程，称为常微分方程；自变量数量2个或以上时，称为偏微分方程。 绝大多数实际工程问题，常微分方程的自变量都是时间t，通常表达为：
F ( t , y , y ˙ , ⋯ , y ( n − 1 ) , y ( n ) ) = 0 Fleft( {t,y,dot y, cdots ,{y^{left( {n - 1} ight)}},{y^{left( n ight)}}} ight) = 0 F(t,y,y˙,⋯,y(n−1),y(n))=0
2) 一阶常微分方程求解上述关于 y ( n ) y^{left( n ight)} y(n)的方程，得到 y ( n ) y^{left( n ight)} y(n)的表达式为：
y ( n ) = f ( t , y , y ˙ , ⋯ , y ( n − 1 ) ) {y^{left( n ight)}} = fleft( {t,y,dot y, cdots ,{y^{left( {n - 1} ight)}}} ight) y(n)=f(t,y,y˙,⋯,y(n−1))
写成向量形式：
[ y ˙ y ¨ ⋮ y ( n ) ] = [ y ¨ y ( 3 ) ⋮ f ( t , y , y ˙ , ⋯ , y ( n − 1 ) ) ] left[ egin{matrix} {dot{y}} \ {ddot{y}} \ vdots \ {{y}^{left( n ight)}} \ end{matrix} ight]=left[ egin{matrix} {ddot{y}} \ {{y}^{left( 3 ight)}} \ vdots \ fleft( t,y,dot{y},cdots ,{{y}^{left( n-1 ight)}} ight) \ end{matrix} ight] ⎣⎢⎢⎢⎡y˙y¨⋮y(n)⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡y¨y(3)⋮f(t,y,y˙,⋯,y(n−1))⎦⎥⎥⎥⎤
令：
y = [ y y ˙ ⋮ y ( n − 1 ) ] mathbf{y}=left[ egin{matrix} y \ {dot{y}} \ vdots \ {{y}^{left( n-1 ight)}} \ end{matrix} ight] y=⎣⎢⎢⎢⎡yy˙⋮y(n−1)⎦⎥⎥⎥⎤
f ( t , y ) = [ y ¨ y ( 3 ) ⋮ f ( t , y , y ˙ , ⋯ , y ( n − 1 ) ) ] = [ y ¨ y ( 3 ) ⋮ f ( t , y ) ] fleft( t,mathbf{y} ight)=left[ egin{matrix} {ddot{y}} \ {{y}^{left( 3 ight)}} \ vdots \ fleft( t,y,dot{y},cdots ,{{y}^{left( n-1 ight)}} ight) \ end{matrix} ight]=left[ egin{matrix} {ddot{y}} \ {{y}^{left( 3 ight)}} \ vdots \ fleft( t,mathbf{y} ight) \ end{matrix} ight] f(t,y)=⎣⎢⎢⎢⎡y¨y(3)⋮f(t,y,y˙,⋯,y(n−1))⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡y¨y(3)⋮f(t,y)⎦⎥⎥⎥⎤
得到向量形式的一阶常微分方程： y ˙ = f ( t , y ) mathbf{dot{y}}=fleft( t,mathbf{y} ight) y˙=f(t,y)
3) 微分方程的解析解法和数值解法微分方程的求解方法有解析解法和数值解法，解析法是求出因变量 y {f{y}} y 关于时间 t t t的具体函数式，表达为 y = y ( t ) {f{y}} = {f{y}}left( t ight) y=y(t)；数值法是解出因变量 y {f{y}} y关于时间 t t t的离散序列，通常表达为 { t ( k ) , y ( k ) } left{ {t(k),{f{y}}(k)} ight} {t(k),y(k)}。 绝大多数的非线性常微分方程，不存在或难以求出解析解，大多数情况下只能求取微分方程的数值解。
2. 算法欧拉法(Euler)是求解一阶常微分方程初值问题的数值方法，分为显式欧拉法、隐式欧拉法、两步欧拉法和改进欧拉法。
2.1 显式欧拉法一阶微分方程初始问题：
{ y ˙ = f ( t , y ) y ( t 0 ) = y 0 left{ egin{array}{l} {f{dot y}} = fleft( {t,{f{y}}} ight) \ {f{y}}left( {{t_0}} ight) = {{f{y}}_0} \ end{array} ight. {y˙=f(t,y)y(t0)=y0
式中， t 0 {t_0} t0为初始时间（已知常数）， y 0 {{f{y}}_0} y0为初始状态（已知向量）， f ( t , y ) fleft( {t,{f{y}}} ight) f(t,y)为关于时间 t {t} t和状态 y {{f{y}}} y的函数（已知函数）。 使用一阶向前差商代替微分，即：
y ˙ ≈ y ( k + 1 ) − y ( k ) t ( k + 1 ) − t ( k ) = y ( k + 1 ) − y ( k ) h {f{dot y}} approx frac{{{f{y}}left( {k + 1} ight) - {f{y}}left( k ight)}}{{tleft( {k + 1} ight) - tleft( k ight)}} = frac{{{f{y}}left( {k + 1} ight) - {f{y}}left( k ight)}}{h} y˙≈t(k+1)−t(k)y(k+1)−y(k)=hy(k+1)−y(k)
式中，h为时间步长。 微分方程变为显式差分方程：
{ y ( k + 1 ) = y ( k ) + h f ( t ( k ) , y ( k ) ) y ( 0 ) = y 0 left{ egin{array}{l} {f{y}}left( {k + 1} ight) = {f{y}}left( k ight) + hfleft( {tleft( k ight),{f{y}}left( k ight)} ight) \ {f{y}}left( 0 ight) = {{f{y}}_0} \ end{array} ight. {y(k+1)=y(k)+hf(t(k),y(k))y(0)=y0
上式是关于 y ( k − 1 ) {f{y}}left( k-1 ight) y(k−1)向 y ( k ) {f{y}}left( k ight) y(k)的递推形式，可以根据初始条件按照递推关系依次求出 y ( 1 ) , y ( 2 ) , y ( 3 ) , y ( 4 ) ⋯ , y ( N ) ⋯ {f{y}}left( 1 ight),{f{y}}left( 2 ight),{f{y}}left( 3 ight),{f{y}}left( 4 ight) cdots ,{f{y}}left( N ight) cdots y(1),y(2),y(3),y(4)⋯,y(N)⋯，此离散序列即为微分方程的数值解。
2.2 隐式欧拉法使用一阶向后差商代替微分，即
y ˙ ≈ y ( k ) − y ( k − 1 ) t ( k ) − t ( k − 1 ) = y ( k ) − y ( k − 1 ) h {f{dot y}} approx frac{{{f{y}}left( k ight) - {f{y}}left( {k - 1} ight)}}{{tleft( k ight) - tleft( {k - 1} ight)}} = frac{{{f{y}}left( k ight) - {f{y}}left( {k - 1} ight)}}{h} y˙≈t(k)−t(k−1)y(k)−y(k−1)=hy(k)−y(k−1)
微分方程变为隐式差分方程：
{ y ( k ) = y ( k − 1 ) + h f ( t ( k ) , y ( k ) ) y ( 0 ) = y 0 left{ egin{array}{l} {f{y}}left( k ight) = {f{y}}left( {k - 1} ight) + hfleft( {tleft( k ight),{f{y}}left( k ight)} ight) \ {f{y}}left( 0 ight) = {{f{y}}_0} \ end{array} ight. {y(k)=y(k−1)+hf(t(k),y(k))y(0)=y0
在第 k k k步迭代时， y ( k − 1 ) , t ( k ) {f{y}}left( {k - 1} ight),t(k) y(k−1),t(k)已知， y ( k ) {f{y}}left( k ight) y(k)为待求未知量， f ( t ( k ) , y ( k ) ) fleft( {tleft( k ight),{f{y}}left( k ight)} ight) f(t(k),y(k))是关于待求未知量 y ( k ) {f{y}}left( k ight) y(k)的函数，即方程为关于 y ( k ) {f{y}}left( k ight) y(k)的非线性方程。通过非线性方程的迭代解法（如牛顿迭代法）可求解出 y ( k ) {f{y}}left( k ight) y(k)。依次求出 y ( 1 ) , y ( 2 ) , y ( 3 ) , y ( 4 ) ⋯ , y ( N ) ⋯ {f{y}}left( 1 ight),{f{y}}left( 2 ight),{f{y}}left( 3 ight),{f{y}}left( 4 ight) cdots ,{f{y}}left( N ight) cdots y(1),y(2),y(3),y(4)⋯,y(N)⋯，此离散序列即为微分方程的数值解。
2.3 两步欧拉法使用中心差商代替微分：
y ˙ ≈ y ( k + 1 ) − y ( k − 1 ) t ( k + 1 ) − t ( k − 1 ) = y ( k + 1 ) − y ( k − 1 ) 2 h {f{dot y}} approx frac{{{f{y}}left( {k + 1} ight) - {f{y}}left( {k - 1} ight)}}{{tleft( {k + 1} ight) - tleft( {k - 1} ight)}} = frac{{{f{y}}left( {k + 1} ight) - {f{y}}left( {k - 1} ight)}}{{2h}} y˙≈t(k+1)−t(k−1)y(k+1)−y(k−1)=2hy(k+1)−y(k−1)
微分方程变为显式差分方程：
{ y ( k + 1 ) = y ( k − 1 ) + 2 h f ( t ( k ) , y ( k ) ) y ( 0 ) = y 0 left{ egin{array}{l} {f{y}}left( {k + 1} ight) = {f{y}}left( {k - 1} ight) + 2hfleft( {tleft( k ight),{f{y}}left( k ight)} ight) \ {f{y}}left( 0 ight) = {{f{y}}_0} \ end{array} ight. {y(k+1)=y(k−1)+2hf(t(k),y(k))y(0)=y0
上式是关于 y ( k − 1 ) {f{y}}left( k-1 ight) y(k−1)， y ( k ) {f{y}}left( k ight) y(k)向 y ( k + 1 ) {f{y}}left( k+1 ight) y(k+1)的递推形式，可以根据初始条件按照递推关系依次求出 y ( 1 ) , y ( 2 ) , y ( 3 ) , y ( 4 ) ⋯ , y ( N ) ⋯ {f{y}}left( 1 ight),{f{y}}left( 2 ight),{f{y}}left( 3 ight),{f{y}}left( 4 ight) cdots ,{f{y}}left( N ight) cdots y(1),y(2),y(3),y(4)⋯,y(N)⋯，此离散序列即为微分方程的数值解。
这里需要特别注意 y ( 1 ) {f{y}}left( 1 ight) y(1)的计算方法：按照公式计算 y ( 1 ) {f{y}}left( 1 ight) y(1)需要用到无意义的 y ( − 1 ) {f{y}}left( -1 ight) y(−1)，此步需要转换成显式欧拉法计算：
{ y ( 1 ) = y ( 0 ) + h f ( t ( 0 ) , y ( 0 ) ) y ( 0 ) = y 0 left{ egin{array}{l} {f{y}}left( 1 ight) = {f{y}}left( 0 ight) + hfleft( {tleft( 0 ight),{f{y}}left( 0 ight)} ight) \ {f{y}}left( 0 ight) = {{f{y}}_0} \ end{array} ight. {y(1)=y(0)+hf(t(0),y(0))y(0)=y0
进一步分析：
y ( k + 1 ) = y ( k − 1 ) + 2 h f ( t ( k ) , y ( k ) ) = [ y ( k − 1 ) + h f ( t ( k ) , y ( k ) ) ] + h f ( t ( k ) , y ( k ) ) {f{y}}left( {k + 1} ight) = {f{y}}left( {k - 1} ight) + 2hfleft( {tleft( k ight),{f{y}}left( k ight)} ight) = left[ {{f{y}}left( {k - 1} ight) + hfleft( {tleft( k ight),{f{y}}left( k ight)} ight)} ight] + hfleft( {tleft( k ight),{f{y}}left( k ight)} ight) y(k+1)=y(k−1)+2hf(t(k),y(k))=[y(k−1)+hf(t(k),y(k))]+hf(t(k),y(k))
式中， y ˉ ( k ) = y ( k − 1 ) + h f ( t ( k ) , y ( k ) ) {f{ar y}}left( k ight) = {f{y}}left( {k - 1} ight) + hfleft( {tleft( k ight),{f{y}}left( k ight)} ight) yˉ(k)=y(k−1)+hf(t(k),y(k)) 相当于隐式欧拉法由 y ( k − 1 ) {f{y}}left( k-1 ight) y(k−1)求解 y ( k ) {f{y}}left( k ight) y(k)； y ( k + 1 ) = y ˉ ( k ) + f ( t ( k ) , y ( k ) ) {f{y}}left( {k + 1} ight) = {f{ar y}}left( k ight) + fleft( {tleft( k ight),{f{y}}left( k ight)} ight) y(k+1)=yˉ(k)+f(t(k),y(k))相当于显式欧拉法由 y ( k ) {f{y}}left( k ight) y(k)求解 y ( k + 1 ) {f{y}}left( k+1 ight) y(k+1)。两步欧拉法是隐式欧拉法与显式欧拉法相结合的一种算法。
2.4 改进欧拉法先使用显式欧拉法求解 k + 1 k+1 k+1步的预测值：
y ˉ ( k + 1 ) = y ( k ) + h f ( t ( k ) , y ( k ) ) {f{ar y}}left( {k + 1} ight) = {f{y}}left( k ight) + hfleft( {tleft( k ight),{f{y}}left( k ight)} ight) yˉ(k+1)=y(k)+hf(t(k),y(k))
再使用 k k k步的微分值 f ( t ( k ) , y ( k ) ) fleft( {tleft( k ight),{f{y}}left( k ight)} ight) f(t(k),y(k))和 k + 1 k+1 k+1步预测的微分值 f ( t ( k + 1 ) , y ˉ ( k + 1 ) ) fleft( {tleft( {k + 1} ight),{f{ar y}}left( {k + 1} ight)} ight) f(t(k+1),yˉ(k+1))的平均值来近似微分：
y ˙ ≈ f ( t ( k ) , y ( k ) ) + f ( t ( k + 1 ) , y ˉ ( k + 1 ) ) 2 {f{dot y}} approx frac{{fleft( {tleft( k ight),{f{y}}left( k ight)} ight) + fleft( {tleft( {k + 1} ight),{f{ar y}}left( {k + 1} ight)} ight)}}{2} y˙≈2f(t(k),y(k))+f(t(k+1),yˉ(k+1))
再利用显式欧拉法校正 k + 1 k+1 k+1步的计算值：
y ( k + 1 ) = y ( k ) + h f ( t ( k ) , y ( k ) ) + f ( t ( k + 1 ) , y ˉ ( k + 1 ) ) 2 {f{y}}left( {k + 1} ight) = {f{y}}left( k ight) + hfrac{{fleft( {tleft( k ight),{f{y}}left( k ight)} ight) + fleft( {tleft( {k + 1} ight),{f{ar y}}left( {k + 1} ight)} ight)}}{2} y(k+1)=y(k)+h2f(t(k),y(k))+f(t(k+1),yˉ(k+1))
综上所述，改进欧拉法的算法为：
{ y ˉ ( k + 1 ) = y ( k ) + h f ( t ( k ) , y ( k ) ) y ( k + 1 ) = y ( k ) + h f ( t ( k ) , y ( k ) ) + f ( t ( k + 1 ) , y ˉ ( k + 1 ) ) 2 y ( 0 ) = y 0 left{ egin{array}{l} {f{ar y}}left( {k + 1} ight) = {f{y}}left( k ight) + hfleft( {tleft( k ight),{f{y}}left( k ight)} ight) \ {f{y}}left( {k + 1} ight) = {f{y}}left( k ight) + hfrac{{fleft( {tleft( k ight),{f{y}}left( k ight)} ight) + fleft( {tleft( {k + 1} ight),{f{ar y}}left( {k + 1} ight)} ight)}}{2} \ {f{y}}left( 0 ight) = {{f{y}}_0} \ end{array} ight. ⎩⎨⎧yˉ(k+1)=y(k)+hf(t(k),y(k))y(k+1)=y(k)+h2f(t(k),y(k))+f(t(k+1),yˉ(k+1))y(0)=y0
上式是关于 y ( k ) {f{y}}left( k ight) y(k)向 y ( k + 1 ) {f{y}}left( k+1 ight) y(k+1)的递推形式，可以根据初始条件按照递推关系依次求出 y ( 1 ) , y ( 2 ) , y ( 3 ) , y ( 4 ) ⋯ , y ( N ) ⋯ {f{y}}left( 1 ight),{f{y}}left( 2 ight),{f{y}}left( 3 ight),{f{y}}left( 4 ight) cdots ,{f{y}}left( N ight) cdots y(1),y(2),y(3),y(4)⋯,y(N)⋯，此离散序列即为微分方程的数值解。
2.5 四种欧拉方法的对比微分方程的数值解法，本质是使用数值积分来实现 y ˙ {f{dot y}} y˙向 y {f{y}} y的转换。显式欧拉法和隐式欧拉法，使用的是矩形数值积分方法，矩形积分是一阶算法，其截断误差为 O ( h 2 ) Oleft( {{h^2}} ight) O(h2)；两步欧拉法和改进欧拉法，使用的是梯形数值积分方法，梯形积分是二阶算法，其截断误差为 O ( h 3 ) Oleft( {{h^3}} ight) O(h3)。
算法类型计算量计算精度稳定性显式欧拉法小低差隐式欧拉法大低差两步欧拉法中较高中改进欧拉法中高好 3. 程序作者使用Matlab开发了四种欧拉方法求解常微分方程的程序，能够方便快捷的求解一阶常微分方程的初值问题。
function [T,X,dX] = ODE_ExplicitEuler( Hfun,t,h,x0 )% [T,X,dX] = ODE_ExplicitEuler( Hfun,t,h,x0 ) 显式欧拉法求解常微分方程% Hfun为描述一阶微分方程导数的函数句柄，格式为 dX = Hfun( t,X )% t为时间节点，可为标量，时间范围为 T = 0:h:t% 长2向量 ，时间范围为 T = t(1):h:t(2)% 向量 ，时间范围为 T = t% h为时间步长，在t的前两种情况下，必须给出h具体值% 直接给出时间节点t时，h可用[]或任意值占位% x0为状态量初始值% 算法：% X(k) = X(k-1) + h*dX(k-1)% By ZFS@wust 2023% 获取