MATLAB使用速成 第八章（定积分的计算与微分方程求解）

一、计算定积分的相关函数

1、trapz（梯形法）

（1）梯形法求解定积分的步骤：

①曲边小梯形的面积可以由直边小梯形的面积来近似，那么整个曲边梯形的面积S为

②将区间[a,b]进行n等分，即令

③整个曲边梯形的面积为

（2）对于语句“trapz(x, y)”，x为分割点（节点）组成的向量，y为被积函数在节点上的函数值组成的向量，它可计算x、y代表的区域的面积，也就是使用梯形法求定积分。另外需要注意的是，其中涉及的运算必须采用数组运算。

2、quad（自适应抛物线法）

（1）自适应抛物线法求解定积分的步骤：

①如下图所示，n等分区间[a,b]，得

②计算节点和中点上的函数值：

③在区间上，用过以下三点的抛物线来近似原函数f(x)。

④设过以上三点的抛物线方程为

⑤在区间上，有

⑥最终求得f(x)在区间[a,b]上的定积分为

（2）对于语句“quad(f,a,b,tol)”，f = f(x)为被积函数（用函数句柄表示），[a,b]为积分区间，tol为计算精度（若不指定，缺省精度是10-6；精度越高，函数运行的时间越长）。另外需要注意的是，其中涉及的运算必须采用数组运算。

3、integral（全局自适应积分法，适用于R2012a以后的版本）

（1）函数integral可用于求解一重积分，具体使用方法为“integral(f,a,b,'RelTol',tol)”，其中f为被积函数的函数句柄（函数涉及的运算必须采用数组运算），[a,b]为积分区间，tol为计算精度（若不指定，缺省精度是10-6，不指定tol时'RelTol'需省略）。

（2）函数integral2可用于求解二重积分，具体使用方法为“integral2(f,a,b,c,d,'RelTol',tol)”，其中f为被积函数的函数句柄（函数涉及的运算必须采用数组运算），[a,b]为第一次积分区间，[c,d]为第二次积分区间，tol为计算精度（若不指定，缺省精度是10-6，不指定tol时'RelTol'需省略）。

二、微分方程求解

（1）对于语句“[T,Y] = solver(odefun,tspan,y0)”：

①y0为初值条件，tspan为求解区间。

②solver为MATLAB的ODE求解器（可以是ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）。

③odefun为显式常微分方程（可理解为一阶导数的函数式），可以用命令inline定义，或在函数文件中定义，然后通过函数句柄调用。

④MATLAB在数值求解时自动对求解区间进行分割，T（列向量）中返回的是分割点的值（自变量），Y（数组）中返回的是这些分割点上的近似解，其列数等于因变量的个数。

（2）没有一种算法可以有效地解决所有的ODE问题，因此MATLAB提供了多种ODE求解器，对于不同的ODE，可以调用不同的求解器。

（3）如果需求解的问题是高阶常微分方程，则需将其化为一阶常微分方程组，此时必须用函数文件来定义该常微分方程组。

①先编写函数文件verderpol.m。

function xprime=verderpol(t,x)
global mu;
xprime=[x(2); mu*(1-x(1)^2)*x(2) - x(1)];

②再编写脚本文件vdpl.m，然后在命令窗口直接运行该文件。

clear;
global mu;
mu=7; y0=[1; 0];
[t,x]=ode45(@verderpol, [0,40], y0);
plot(t,x(:,1));