第7次课 数值微积分与常微分方程求解-修订

6.3.1 龙格-库塔法简介
对于一阶常微分方程的初值问题，在求解未知函数 y 时，y 在 t0 点的值 y(t0 ) = y0 是已知的， 并且根据高等数学中的中值定理，应有 y(t0 + h) = y1 ?y0 + hf(t0 ,y0 )，h>0，称为步长 y(t0 + 2h) = y2 ?y1 + hf(t1 ,y1 ) 一般地，在任意点 ti = t0 +ih，有 y(t0 + ih) = yi ?yi? 1 + hf(ti? 1 ,yi?1 )，i = 1，2，„，n 当然，递推过程中有一个误差累计的问题。在实际计算过程中，使用的递推公式一般进行过 改造，著名的 4 阶龙格-库塔公式为
表 6.1 求 解 器
求常微分方程数值解的函数 采 用 方 法 2-3 阶 Runge-Kutta 算法，低精度。 4-5 阶 Runge-Kutta 算法，中精度。 Adams 算法，精度可到 10 -3 ～10 -6 梯形算法 Gear’s 反向数值微分算法，中精度。 2 阶 Rosebrock 算法，低精度。 梯形算法，低精度。 可变秩方法 非刚性 非刚性 非刚性，计算时间比 ode45 短 适度刚性 刚性 刚性，当精度较低时，计算时间比 ode15s 短 刚性，当精度较低时，计算时间比 ode15s 短 完全隐式微分方程 适 用 场 合
一般情况下，quadl函数调用的步数明显小于quad函数，而且精度更高， 1 从而保证能以更高的效率求出所需的定积分值。 1 dx 【例】分别用quad函数和quadl函数求椭圆积分 0 1 x 4 的近似值， 并在相同的积分精度下，比较函数的调用次数。 调用函数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
gx=cos(x); plot(x,dpx,x,dx,'o',x,gx,'+');
定积分的数值求解实现
1．自适应辛普生法 quad函数和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为被积函数 的调用次数。
【例】设f(x) = sinx，用不同的方法求函数f(x)的数值导数， 并在同一个坐标系中做出f '(x)的图像。 为确定计算数值导数的点，假设在[0， π]区间内以π/24为步 长求数值导数。下面用3种方法求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); dx=diff(sin([x,pi+pi/24]))/(pi/24);
称f(x)、f(x)及 δf(x)分别为函数在 x 点处以 h（h >0）为步长的向前差分、向后差分和中心差分。 当步长 h 充分小时，有
f '( x)
f ( x ) h
f '( x) f '( x)
f ( x ) h
f ( x)
h
和差分一样，称f(x)/h、f(x)/h 及 δf(x)/h 分别为函数在 x 点处以 h（h>0） 为步长的向前差商、 向后差商和中心差商。当步长 h（h >0）充分小时，函数 f 在点 x 的微分接近于函数在该点的任意 种差分，而 f 在点 x 的导数接近于函数在该点的任意种差商。
ode23 ode45 ode113 ode23t ode15s ode23s ode23tb ode15i
若微分方程描述的一个变化过程包含着多个相互作用但变化速度相差十 分悬殊的子过程，这样一类过程就认为具有“刚性”，这类方程具有非 常分散的特征值。求解刚性方程的初值问题的解析解是困难的，常采用 表6.1中的函数ode15s、ode23s和ode23tb求其数值解。求解非刚性的一 阶常微分方程或方程组的初值问题的数值解常采用函数ode23和ode45， 其中ode23采用2阶龙格-库塔算法，用3阶公式作误差估计来调节步长， 具有低等的精度；ode45则采用4阶龙格-库塔算法，用5阶公式作误差估 计来调节步长，具有中等的精度。
k1 f ( xi , yi ) k2 f ( x 1 , yi h k1 ) i 2 2 h k3 f ( xi 1 , yi k2 ) 2 2 k f ( x , y hk ) i 1 i 3 4 h yi 1 yi (k1 2k2 2k3 k4 ) 6
2
常微分方程的数值求解
科学实验和生产实践中的许多原理和规律都可以描述成适当 的常微分方程，如牛顿运动定律、生态种群竞争、股票的 涨幅趋势、市场均衡价格的变化等，许多情况下都要进行 数值求解。 常微分方程初值问题的数值解法，首先是建立求数值解的递 推公式。递推公式通常有两类，一类是计算yi +1时只用到 xi +1、xi和yi，即前一步的值，此类方法称为单步法，其 代表是龙格-库塔（Runge-Kutta）法。另一类是计算yi+1 时，需要前面k步的值，此类方法称为多步法，其代表是 亚当姆斯（Adams）法。这些方法都是基于把一个连续的 定解问题离散化为一个差分方程来求解，是一种步进式的 方法。
e 法 函数trapz 该函数调用格式如下。 ● T = trapz(Y)：若Y是一向量，则从1开始取单位步长，以 Y的值为函数值计算积分值。若Y是一矩阵，则计算Y的每 一列的积分。 ● T = trapz(X,Y)：向量X、Y定义函数关系Y = f(X)。X、Y 是两个等长的向量：X = (x1，x2，…，xn)，Y = (y1， y2，…，yn)，并且x1<x2<…<xn，积分区间是[x1，xn]。
【例】生成一个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) 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列 是等间距的。
options用命令odeset生成设置求解属性常用的属性包括相对误差值?reltol?默认时为10与绝对误差向量?abstol?默认时每一元素为10若微分方程描述的一个变化过程包含着多个相互作用但变化速度相差十分悬殊的子过程这样一类过程就认为具有刚性这类方程具有非常分散的特征值
数值微积分与常微分方程求解
多重定积分的数值求解实现 多重积分的被积函数是二元函数或三元函数，积分范围是平 面上的一个区域或空间中的一个区域。MATLAB中提供的 d b dblquad函数用于求 的数值解， c a f ( x, y )dxdy triplequad函数用于求 f d b f ( x, y, z )dxdydz 的数值解。函数 e c a 的调用格式为 ● 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完全相同。
【例】计算二重定积分 1 2 e x / 2 sin( x 2 y )dxdy 1 2 （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)
数值微分的实现 1．基本思想是先用逼近或拟合等方法将已知数据在一定范围内的近似函 数求出，再对此近似函数进行微分。 比如用多项式函数g(x)对f(x)进行逼近，然后用g(x)在点x处的导数作 为f(x)在点x处的导数。该种方法一般只用在低阶数值微分。 2．用f(x)在点x处的某种差商作为其导数。 在MATLAB中，提供计算向前差分的函数diff。 ● DX = diff(X)：计算向量X的向前差分，DX(i) = X(i+1) −X(i)， i = 1,2, …,n −1。 ● DX = diff(X,n)：计算X的n阶向前差分。例如，diff(X,2) = diff(diff(X))。 ● DX = diff(A,n,dim)：计算矩阵A的n阶差分，dim = 1时（默认状态）， 按列计算差分；dim = 2，按行计算差分。 函数diff计算的是向量元素间的差分，故所得输出比原向量少了一个 元素。
2．高斯-克朗罗德法
该函数的调用格式为 [I,err] = quadgk(@fname,a,b) 其中，err返回近似误差范围，其他参数的含义和用法与quad函数相同。 积分上下限可以是−Inf或Inf，也可以是复数。如果积分上下限是复数， 则quadgk在复平面上求积分。 【例】求定积分 0 （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个括号里面是代表函数运算的表达式。
【学习目标】 ● 掌握微分与积分的数值计算方法。 ● 掌握常微分方程的数值求解方法。
数值微分
数值差分与差商
f ( x) f ( x h) f ( x) f ( x ) f ( x ) f ( x h )
f ( x) f ( x h / 2) f ( x h / 2)
4 阶标准龙格-库塔法精度高，程序简单，易于改变步长，比较稳定，是一个常用的方法，但 计算量较大。当函数 f(x，y)较为复杂，可用显式亚当姆斯方法或亚当姆斯预测 —校正方法，不仅 计算量较小，稳定性也比较好，但不易改变步长。
常微分方程数值求解的实现 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'（默认时为10−3）与绝对误差 向量'AbsTol'（默认时每一元素为10−6）。