第七章微积分的数值计算

1、7.1 数值微分7.2 数值积分7.3 常微分方程的数值解法 7.1 数值微分数值微分 实际问题常需计算函数的导数或积分值。但很多情况下，函数关系难以准确表示；即使能使用解析式准确表示，表示式却很复杂，不能用于实际计算。本章介绍数值计算导数或积分的实用方法。7.1.1 差分和差商差分和差商 根据导数的定义其中，x和y分别称为自变量x和因变量y的增量，也称之为差分。可以用差分的商 作微商（导数）的近似。数值微分就是用函数值的线性组合近似函数在某点的导数值。自变量x的步长一般取定值。 首先在xi处对函数进行泰勒展开，1001limlimiiixxiiyyyyxxx yx234( 4 )12 !3

2、!4 !iiiiiixxxyyxyyyy 根据不同的组合方式可以得到精度不同的差分公式。以函数的一阶导数为例 ：微分公式微分公式表达式表达式截断误差截断误差两点前向o(x)两点后向o(x)三点中心o(x2)三点前向o(x2)三点后向o(x2)五点中心o(x4)1iiiyyyx 1iiiyyyx 112iiiyyyx 12432iiiiyyyyx12432iiiiyyyyx 11228812iiiiiyyyyyx精度为o(x2)的高阶中心差分算法11221123(4)211242222464iiiiiiiiiiiiiiiyyyyxyyyyyxyyyyyyx精度为o(x4)的高阶中心差分算法211

3、2211223211233(4)3211234881216301612813138812395639126iiiiiiiiiiiiiiiiiiiiiiiiiiyyyyyxyyyyyyxyyyyyyyxyyyyyyyyx7.1.2 数值微分的实现数值微分的实现 在matlab中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff和梯度函数gradient。diff调用格式为：调用格式为： dy=diff(y)：计算向量y的向前差分，并把结果赋值给向量dy dy(i)=y(i+1)-y(i)，i=1,2,n-1。注意向量dy元素个数比y少dy=diff(y,n)：计算向量y的n阶向前差分

4、。例如， diff(y,2)=diff(diff(y)=dx(i+1)-dx(i)= y(i+2)-2y(i+1)+y(i) ， i=1,2 n-2。dx=diff(a,n,dim)：计算矩阵a的n阶差分，dim=1时(缺省状态)，按列计算差分；dim=2，按行计算差分。 a=pascal(4)a = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 b=diff(a)b = 0 1 2 3 0 1 3 6 0 1 4 10 c=diff(a,2)c = 0 0 1 3 0 0 1 4 d=diff(a,1,1)d = 0 1 2 3 0 1 3 6 0 1 4 10 e=

5、diff(a,1,2)e = 0 0 0 1 1 1 2 3 4 3 6 107.1.3 近似梯度函数近似梯度函数gradient的调用格式为的调用格式为 fx,fy=gradient(f,hx,hy) 求矩阵f， 求其x(行)方向的数值梯度fx和y(列)方向的数值梯度fy，x方向步长全为hx，y方向步长全为hy。 fx相当于偏导数f/x ,fy相当于偏导数f/y。 fx,fy=gradient(f,h) 求矩阵f， 求其x(行)方向的数值梯度fx和y(列)方向的数值梯度fy，各个方向步长全为h。 fx=gradient(f) 如果f是向量，直接求其数值梯度；如果f是矩阵，求其x(行)方向的数

6、值梯度，步长为1 。 fx,fy=gradient(f) 求矩阵f， 求其x(行)方向的数值梯度fx和y(列)方向的数值梯度fy，步长为1ffffijkxyz x=1 3 5 2 4;y= gradient(x)y = 2.0000 2.0000 -0.5000 -0.5000 2.0000 y=gradient(x,2)y= 1.0000 1.0000 -0.2500 -0.2500 1.0000即两两边用前向和后向差分，中间用中边用前向和后向差分，中间用中心差分心差分 a=pascal(4)a = 1 1 1 1 1 2 3 4 1 3 6 10 1 4 10 20 fx,fy=gradi

7、ent(a)fx = 0 0 0 0 1.0000 1.0000 1.0000 1.0000 2.0000 2.5000 3.5000 4.0000 3.0000 4.5000 8.0000 10.0000fy = 0 1.0000 2.0000 3.0000 0 1.0000 2.5000 4.5000 0 1.0000 3.5000 8.0000 0 1.0000 4.0000 10.0000211315342234545,222xxyhxxxxxxxyyyhxhxhxxxyhx，7.1.4 拉普拉斯算子拉普拉斯算子4*del2 由于内部算法的原因： u=4*del2(v,h)，对1维向量

8、v以步长h求拉普拉斯算子时，返回一相同维数的向量u，且默认的步长为1。 u=4*del2(v,h1,h2)，对矩阵v，横向(x方向)以步长h1，纵向(y方向)以步长 h2计算拉普拉斯算子。2222211( ,)( ,)( ,)44vx yvx yuvx yxy14321211212321(452)1(2)1( 254)iiiinnnnnuvvvvhuvvvhuvvvvh 4*del2(u)ans = 4 4 4 4 4 4 4 4 4 4 4 4x,y= meshgrid(1:4,1:3); u = x.*x+y.*yu = 2 5 10 17 5 8 13 20 10 13 18 257.2

9、 数值积分数值积分7.2.1 数值积分基本原理数值积分基本原理 我 们 知 道 ， 定 积 分 是 求 和 式 的 极 限 ，即 。它的几何意义是曲边梯形的面积。从定义可知，定积分的基本分析方法是四步，即分割、近似、求和、取极限。分割就是把总量（整块曲边梯形面积）分成若干分量（小曲边梯形面积）；近似就是在每个分量中用容易计算的量去代表；求和就是把分量加起来得到总近似值；最后取极限就得到积分精确值。 1( )dlim()nbkkankfxxfx把区间a，b分割成n等分，像这样取定步长算积分的方法，称为定步长积分法法。常见的低阶求积分公式常见的低阶求积分公式 复化矩形公式 复化的梯形公式 复化的辛

10、普森（simpson）公式 10( )bnkkaf x dxs1/21()4 ()()6kkkkhsf xf xf x()kksfxh1()()2kkkhsfxfx, (0,1,1)kbaxakh hknn辛普森公式的几何意义辛普森公式的几何意义 7.2.2 变步长积分法变步长积分法 计算积分，可以采取逐步缩小步长h的办法。即先任取步长h进行计算，然后取较小步长 h 进行计算，如果两次计算结果相差较大，则取更小步长进行计算，如此下去，直到相邻两次计算结果相差不大为止，取最小步长算出的结果作为积分值。这种方法称为变步长积分法。 利用两种步长计算积分时，通常取h=h/2 。而每次改变步长后，只需计

11、算新增节点处的函数值，将它们的和乘新步长 。 matlab常用的数值积分函数参数名参数名功能描述功能描述quad一元函数的数值积分，采用变步长辛普森simpson法(低阶)quadl一元函数的数值积分，采用变步长洛巴托lobatto法(高阶)qualv一元函数的矢量数值积分dbquad二重积分（默认值为simpson法）triplequad三重积分（默认值为simpson法） 7.2.3 一元函数的数值积分举例一元函数的数值积分举例求定积分1. 建立匿名函数建立匿名函数 f=(x)x./(x.2+1);2. 用用sum函数实现复化矩形法求积函数实现复化矩形法求积1122001lo g (1)0

12、 .3 4 6 612exd xxxx=0:0.001:1; %步长提高到0.001 y=f(x); sum(y)*0.001ans = 0.3468 x=0:0.01:1; %步长0.01 y=f(x); sum(y)*0.01ans = 0.3491 3. 用用trapz函数实现复化梯形法求积函数实现复化梯形法求积x=0:0.001:1;y=f(x);y=trapz(y)*0.001y = 0.3466x=0:0.01:1;y=f(x);y=trapz(y)*0.01y = 0.34664. 用用matlab函数求定积分函数求定积分用quad实现变步长辛普森法求定积分。调用格式为 quad

13、(fun,a,b,tol)其中fun为积分函数，a,b为积分区间，tol为积分的误差阈值，默认值为1e-6quad(f,0,1)ans = 0.3466用quadl实现变步长洛巴托求定积分。调用格式为 quadl(fun,a,b,tol)其中fun为积分函数，a,b为积分区间，tol为积分的误差阈值，默认值为1e-6 quadl(f,0,1)ans = 0.34667.2.4 矢量积分矢量积分相当于多个一元定积分。例如求y=(x,n)1./(sqrt(2*pi).*(1:n).*exp(-x.2./(2*(1:n).2); %归一化高斯函数quadv(x)y(x,5),-1,1)ans = 0

14、.6827 0.3829 0.2611 0.1974 0.1585矢量积分的结果是一个向量，每一个元素为一个一元函数定积分的值。12211exp(),1,2,3,4,522xdx nn7.2.5 二重定积分的数值求解二重定积分的数值求解 使用matlab提供的dblquad函数就可以直接求出上述二重定积分的数值解。该函数的调用格式为： dblquad(f,a,b,c,d,tol,method)该函数求f(x,y)在a,bc,d区域上的二重定积分。参数tol，可以采用method=quadl的方法使函数用高阶的洛巴托法计算定积分。例如计算 f=(x,y)exp(-x.2-y.2)/pi; %归一

15、化高斯函数dblquad(f,-1,1,-1,1,1e-6,quadl)ans = 0.7101三重积分函数triplequad用法与二重积分类似1122111exp()xyd xdy 7.3 常微分方程的数值解法常微分方程的数值解法7.3.1 常微分方程初值问题的数值解法常微分方程初值问题的数值解法一般形式为所谓的数值解法就是求解y(x)在 区间的近似值yn的方法，yn(n=1,2,n)称为常微分方程的数值解。自变量x的步长一般为定值h。0( ,)( )dyfx yaxbdxy ay01na xxxb7.3.2 常见的数值方法常见的数值方法向前欧拉公式向后欧拉公式梯形公式改进的欧拉公式10(

16、,)0,1( )nnnnyyhfxynnyy a1110(,)0,1( )nnnnyyhfxynnyy a1110(,)(,)0,12( )nnnnnnhyyf xyf xynnyy a1(,)(,)1()2pnnnqnnpnpqyyhfxyyyhfxh yyyy 7.3.3 龙格库塔法简介龙格库塔法简介 基本思想就是利用在某些点处的值的线性组合构造公式，使其按泰勒展开后与初值问题的解的泰勒展开式比较，有尽可能多的相同项，从而保证算式有较高的精度。常用的四阶经典的龙格-库塔公式112341213243(22)6(,)(,)22(,)22(,)22nnnnnnnnnnhyykkkkkf x yh

17、hkf xykhhkf xykhhkf xyk7.3.4 龙格库塔法的实现龙格库塔法的实现求解器求解器ode类型类型特点特点精度精度说明说明ode45非刚性一步算法（只需前一步的结果），4,5阶 runge-kutta方法。中大部分场合的首选算法ode23非刚性一步算法，2,3阶runge-kutta方法。低使用于精度较低的情形ode113非刚性多步法（需要前几布的结果），adams算法。低高计算时间比ode45短ode23t适度刚性 采用梯形算法适度刚性情形ode15s刚性多步法，gears反向数值积分。低中若ode45失效时，可尝试使用ode23s刚性一步法，2阶rosebrock算法精度

18、。低当精度较低时，计算时间比ode15s短 对于一个常微分方程组，如果其解相差十分悬殊，就称之为刚性方程组刚性方程组。对于刚性方程组，为了保持解法的稳定，步长选取十分困难，有些解法不再适用。ode函数调用格式：t,y = ode ij (odefun,tspan,y0)t,y = ode ij (odefun,tspan,y0,options)t,y = ode ij (odefun,tspan,y0,options,p1,p2,.)t,y,te,ye,ie=odeij(odefun,tspan,y0,options,p1,p2,.)odefun为显式常微分方程中的 f(xn,yn)，tspan为求解区间，y0为初始条件。7.3.5 求解多变量一阶常微分方程组。求解多变量一阶常微分方程组。例：求解下面的混沌理论的洛仑兹方程。dx8= -y + yzdt3dy= 10y + 10zdtdz= - xy + 35y- zdt编写lorfun.m如下：function ydot=lorfun(t,y)ydot=-8/3*y(1)+y(2)*y