第6章数值微积分与常微分方程求解,在许多实际问题中要采用数值方法来求函数的微分或积分。 而在实际问题中遇到的常微分方程往往很复杂，在许多情况下得不出一般解，所以，一般是要求获得解在若干个点上的近似值。 【本章学习目标】 掌握微分与积分的数值计算方法。 掌握常微分方程的数值求解方法。,6.1 数 值 微 分,6.1.1 数值差分与差商,【例6.1】设f(x)=sinx，用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做出f(x)的图像。 为确定计算数值导数的点，假设在0，20区间内以为步长求数值导数。下面用3种方法求f(x)在这些点的导数。首先用一个5次多项式p(x)拟合函数f(x)，并对p(x)求一般意义下的导数dp(x)，求出dp(x)在假设点的值；第2种方法用diff函数直接求f(x)在假设点的数值导数；第3种方法先求出导函数f(x)=cosx，然后直接求f(x)在假设点的导数。 x=0:pi/24:pi; %用5次多项式p拟合f(x),并对拟合多项式p求导数dp在假设点的函数值 p=polyfit(x,sin(x),5); dp=polyder(p); dpx=polyval(dp,x); %直接对sin(x)求数值导数 dx=diff(sin(x,pi+pi/24)/(pi/24); %求函数f的导函数g在假设点的导数 gx=cos(x); %作图 plot(x,dpx,x,dx,o,x,gx,+);,对于求矩阵的差分，即为求各列或各行向量的差分，从向量的差分值可以判断列或行向量的单调性、是否等间距以及是否有重复的元素。 【例6.2】生成一个5阶魔方矩阵，按列进行差分运算。 M=magic(5) M= 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9 DM=diff(M) %计算V的一阶差分 DM= 6 19 6 6 1 19 1 6 6 6 6 6 6 1 19 1 6 6 19 6 可以看出，diff函数对矩阵的每一列都进行差分运算，因而结果矩阵的列数是不变的，只有行数减1。矩阵DM第3列值相同，表明原矩阵第3列是等间距的。,6.2.2 定积分的数值求解实现 在MATLAB中可以使用quad或quadl来进行数值积分。 1自适应辛普生法 MATLAB提供了基于自适应Simpson法的quad函数和自适应Lobatto法的quadl函数来求定积分。函数的调用格式为 I,n=quad(fname,a,b,tol,trace) I,n=quadl(fname,a,b,tol,trace) 其中fname是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，默认时取tol=10-6。trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认时取trace=0。返回参数I即定积分值，n为被积函数的调用次数。 【例6.3】求S=。 （1）建立被积函数文件fe.m。 function f=fe(x) f=x.*sin(x)./(1+abs(cos(x); （2）调用数值积分函数quad求定积分。 quad(fe,0,pi),一般情况下，quadl函数调用的步数明显小于quad函数，而且精度更高，从而保证能以更高的效率求出所需的定积分值。 【例6.4】分别用quad函数和quadl函数求椭圆积分 的近似值，并在相同的积分精度下，比较函数的调用次数。 调用函数quad求定积分： format long; fx=inline(1./sqrt(1+X.4); %定义一个语句函数 I,n=quad(fx,0,1,1e-10)%注意函数名不加号 I= 0.927037338654481 n= 113 调用函数quadl求定积分： I,n=quadl(fx,0,1,1e-10) I= 0.927037338650659 n= 48,2高斯-克朗罗德法 MATLAB提供了基于自适应高斯-克朗罗德法的quadgk函数来求振荡函数的定积分。该函数的调用格式为 I,err=quadgk(fname,a,b) 其中，err返回近似误差范围，其他参数的含义和用法与quad函数相同。积分上下限可以是Inf或Inf，也可以是复数。如果积分上下限是复数，则quadgk在复平面上求积分。 【例6.5】求定积分 （1）建立被积函数文件feln.m。 function y=feln(x) y=exp(x).*log(x); （2）调用数值积分函数quadgk求定积分。 format long; I=quadgk(feln,0,1) I= 1.317902162414081 也可使用匿名函数方法求解： I=quadgk(x)(exp(x).*log(x),0,1); 匿名函数是构造简单函数的一种方法，是函数指针，第1个括号里面是自变量，第2个括号里面是代表函数运算的表达式。,3梯形积分法 在MATLAB中，对由表格形式定义的函数关系的求定积分问题用梯形积分函数trapz。该函数调用格式如下。 T = trapz(Y)：若Y是一向量，则从1开始取单位步长，以Y的值为函数值计算积分值。若Y是一矩阵，则计算Y的每一列的积分。例如： trapz(1:5;2:6) ans= 12 16 T=trapz(X,Y)：向量X、Y定义函数关系Y=f(X)。X、Y是两个等长的向量：X=(x1，x2，xn)，Y=(y1，y2，yn)，并且x1x2xn，积分区间是x1，xn。,6.2.3 多重定积分的数值求解实现 定积分的被积函数是一元函数，积分范围是一个区间；而重积分的被积函数是二元函数或三元函数，积分范围是平面上的一个区域或空间中的一个区域。MATLAB中提供的dblquad函数用于求 的数值解，triplequad函数用于求 的数值解。函数的调用格式为 dblquad(fun,a,b,c,d,tol,trace) triplequad(fun,a,b,c,d,e,f,tol,trace) 其中，fun为被积函数，a，b为x的积分区域，c，d为y的积分区域，e，f为z的积分区域，参数tol、trace的用法与函数quad完全相同。,【例6.7】计算二重定积分 （1）建立一个函数文件fxy.m： function f=fxy(x,y) global ki; ki=ki+1; %ki用于统计被积函数的调用次数 f=exp(-x.2/2).*sin(x.2+y); （2）调用dblquad函数求解： global ki; ki=0; I=dblquad(fxy,-2,2,-1,1) Ki 如果使用inline函数，则命令如下： f=inline(exp(-x.2/2).*sin(x.2+y),x,y); I=dblquad(f,2,2,1,1),6.3 常微分方程的数值求解,科学实验和生产实践中的许多原理和规律都可以描述成适当的常微分方程，如牛顿运动定律、生态种群竞争、股票的涨幅趋势、市场均衡价格的变化等。 常微分方程初值问题的数值解法，首先要解决是建立求数值解的递推公式。递推公式通常有两类，一类是计算yi+1时只用到xi+1、xi和yi，即前一步的值，此类方法称为单步法，其代表是龙格-库塔（Runge-Kutta）法。另一类是计算yi+1时，需要前面k步的值，此类方法称为多步法，其代表是亚当姆斯（Adams）法。这些方法都是基于把一个连续的定解问题离散化为一个差分方程来求解，是一种步进式的方法。,6.3.1 龙格-库塔法简介,6.3.2 常微分方程数值求解的实现 MATLAB提供了多个求常微分方程数值解的函数，一般调用格式为 t,y=solver(fname,tspan,y0,options) 其中t和y分别给出时间向量和相应的状态向量。solver为求常微分方程数值解的函数ode23、ode45、ode113、ode23t、ode15s、ode23s、ode23tb、ode15i之一，表6.1所示为各函数的采用方法和适用问题。fname是定义f(t，y)的函数文件名，该函数文件必须返回一个列向量。tspan形式为t0，tf，表示求解区间。y0是初始状态列向量。options（用命令odeset生成）设置求解属性，常用的属性包括相对误差值RelTol（默认时为103）与绝对误差向量AbsTol（默认时每一元素为106）。,若微分方程描述的一个变化过程包含着多个相互作用但变化速度相差十分悬殊的子过程，这样一类过程就认为具有“刚性”，这类方程具有非常分散的特征值。求解刚性方程的初值问题的解析解是困难的，常采用表6.1中的函数ode15s、ode23s和ode23tb求其数值解。求解非刚性的一阶常微分方程或方程组的初值问题的数值解常采用函数ode23和ode45，其中ode23采用2阶龙格-库塔算法，用3阶公式作误差估计来调节步长，具有低等的精度；ode45则采用4阶龙格-库塔算法，用5阶公式作误差估计来调节步长，具有中等的精度。,（1）建立函数文件funst.m。 function yp=funst(x,y) yp=sec(x)+y*tan(x); （2）求解微分方程。 x0=0; xf=1; y0=pi/2; x,y=ode23(funst,x0,xf,y0);%求数值解 y1=(x+pi/2)./cos(x); %求精确解 x,y,y1,（1）建立函数文件verderpol.m。 function xprime=verderpol(t,x) global mu; xprime=x(2);mu*(1x(1)2)*x(2)x(1); （2）求解微分方程。 global mu; mu=2; y0=1;0; t,x=ode45(verderpol,0,20,y0); （3）绘制解的曲线。 subplot(1,2,1);plot(t,x); %系统时间响应曲线 subplot(1,2,2);plot(x(:,1),x(:,2) %系统相平面曲线,（1）建立模型的函数文件lorenz.m。 function xdot=lorenz(t,x) xdot=8/3,0,x(2);0,10,10;x(2),28,1*x; （2）解微分方程组。 t,x=ode23(lorenz,0,80,0;0;eps); （3）绘制系统相平面图 plot3(x(:,1),x(:,2),x(:,3); axis(10,45,-15,20,-30,25);,