常微分方程的二阶泰勒数值解法及Matlab程序及案例

扬辰鑫

已于 2024-05-31 14:07:46 修改

阅读量142

点赞数 1

文章标签： matlab 人工智能算法

于 2024-05-16 10:30:00 首次发布

原文链接：https://mp.weixin.qq.com/s/soj6Ce0nLqYMZKej5cEVQw

版权

1.算法

微分方程数值解法的核心是使用数值积分来实现状态导数dy/dt向状态y的转换，欧拉法只利用了状态的一阶导数信息，二阶泰勒方法通过引入二阶导数信息来提高求解精度。

对于一阶微分方程初值问题：

式中，t0为初始时间（已知常数），y0为初始状态（已知向量），f(t,y)为关于时间t和状态y的函数（已知函数）。

将微分方程解y(t)在工作点附近展开了泰勒级数，并忽略二阶以上的高阶项：

由于：

上式可进一步简化为：

定义步长：

按照上式求出：

时的函数值:：

其中导数可以通过对状态y和时间t分别求偏导方式来简化计算：

，则二阶泰勒数值解法为：

在第k步迭代时，y(k)、t(k)已知，y(k+1)为待求未知量，可以根据初始条件按照递推关系依次求出y(1)、y(2)、y(3)…y(k)…，此离散序列即为微分方程的数值解。

2 原理

微分方程的数值解法，本质是使用数值积分来实现dy/dt向y的转换。四阶龙格库塔法通过对微分的四步分段逼近，在一个求解步长内能够逼近复杂的曲线，因此能够取得较高的计算精度，其截断误差为O(h^5)。

注意到递推公式中用到了二阶导数信息，也正是因为利用了二阶导数信息，在一个步长内使用梯形数值积分，相对于欧拉法，在一个求解步长内能够逼近复杂的曲线，其泰勒展式中的截断误差相对欧拉方法更小，其截断误差为O(h^3)，在相同的求解步长条件下，二阶泰勒方法精度会优于欧拉法。二阶泰勒方法的数值积分过程如下图所示。

图1 数值积分示意图

3 程序

作者使用Matlab开发了二阶泰勒法求解常微分方程的程序，能够方便快捷的求解一阶常微分方程的初值问题。

function [T,X,dX] = ODE_Taylor2( Hfun,t,h,x0 )
% [T,X] = ODE_Taylor2( Hfun,t,h,x0 ) 2阶泰勒方法求解常微分方程
% Hfun为描述状态导数的函数句柄，格式为 [dX,d2X] = Hfun( t,X )
% dX为一阶导数，d2X为二阶导数
% 必须保证返回dX的格式为行向量
% 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 + d2X*(h/2) )