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  常微分方程数值