4.8  微分方程

微分方程是数值计算中常见的问题，MATLAB提供了多种函数来计算微分方程的解。

4.8.1  常微分方程

众所周知，对一些典型的常微分方程，能求解出它们的一般表达式，并用初始条件确定表达式中的任意常数。但实际中存在有这种解析解的常微分方程的范围十分狭窄，往往只局限在线性常系数微分方程(含方程组)，以及少数的线性变系数方程。对于更加广泛的、非线性的一般的常微分方程，通常不存在初等函数解析解。由于实际问题求解的需要，求近似的数值解成为了解决问题的主要手段。常见的求数值解的方法有欧拉折线法、阿当姆斯法、龙格-库塔法与吉尔法等。其中由于龙格-库塔法的精度较高，计算量适中，所以使用的较广泛。

数值解的最大优点是不受方程类型的限制，即可以求任何形式常微分方程的特解(在解存在的情况下)，但是求出的解只能是数值解。

1．龙格-库塔方法简介

对于一阶常微分方程的初值问题，在求解未知函数时，在点的值是已知的，并且根据高等数学中的中值定理，应有：

一般而言，在任意点，有：

当确定后，根据上述递推公式能算出未知函数在点的一列数值解：

当然在递推的过程中同样存在着一个误差累计的问题，实际计算中的递推公式一般都进行过改造，龙格-库塔公式为：

2．龙格-库塔法的实现

基于龙格-库塔法，MATLAB提供了ode系列函数求常微分方程的数值解。常用的有ode23 和ode45函数，其调用语法如下。

(1)[t,y]=ode23(filename,tspan,y0)：采用了二阶、三阶龙格-库塔法进行计算。

(2)[t,y]=ode45(filename,tspan,y0)：采用了四阶、五阶龙格-库塔法进行计算。

其中filename是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。tspan的形式为[t0,tf]，表示求解区间。y0是初始状态列向量。t和y分别给出时间向量和相应的状态向量。

这两个函数分别采用了二阶、三阶龙格-库塔法和四阶、五阶龙格-库塔法，并采用自适应变步长的求解方法，即当解的变化较慢时采用较大的步长，从而使得计算的速度很快；当解的变化较快时步长会自动地变小，从而使得计算的精度很高。

【例4-45】  设有初值问题：

试求其数值解，并和精确解相比较，精确解为()。

首先要建立微分方程所对应的函数文件myodefun.m，文件内容如下：

function y=myodefun(t,y)

%  建立函数文件myodefun.m

y=(y^2-t-2)/(4*(t+1));

建立myodefun函数之后，就可以调用ode23函数求解微分方程。

>> t0=0;

>> tf=10;

>> y0=2;

>> [t,y]=ode23 ('myodefun',[t0,tf],y0);       %  求数值解

>> y1=sqrt(t+1)+1;                              %  求精确解

>> plot(t,y,'k.',t,y1,'r')

通过图形比较，数值解用黑色圆点表示，精确解用红色实线表示，如图4-13所示。

【例4-46】  求下面无劲度系统微分方程组的数值解。

为了求解方程，首先要建立方程的m文件。本例中不妨建立名为rigid.m的函数文件，此文件用以描述给出的方程组，文件的内容如下：

function dy = rigid(t,y)

dy = zeros(3,1);                 %  一个列向量

dy(1) = y(2) * y(3);

dy(2) = -y(1) * y(3);

dy(3) = -0.51 * y(1) * y(2);

本例中，我们通过odeset函数对误差进行控制，另外在时间[0 12]进行求解，0时刻初始条件向量为[0 1 1]。

>> options =odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);%  误差控制

>> [T,Y] = ode45(@rigid,[0 12],[0 11],options);               %  求数值解

>>plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')               %  绘制结果图

得到的结果如图4-14所示。

图4-13  常微分方程结果图

图4-14  常微分方程数值