94 第六章 计算方法简介§1 数值逼近1.1插值许多实际问题都要用函数)(x f y =来表示某种内在规律的数量关系，其中相当一部分函数虽然可能在某个区间上具有很好的性质（连续、光滑等），但没有函数的表达式信息，我们只能通过实验或者观测得到函数在一些点i x 上的函数值)(i i x f y =),2,1,0(n i =，这是一张函数表.有些函数虽然有解析式，但由于计算复杂，使用不方便，我们通常也造一个函数表，例如三角函数表、平方根表等.为了研究函数的性质，往往还需要求出不在函数表上的函数值，因此我们希望根据给定的函数表构造一个既能反映函数)(x f y =的性质、又便于计算的简单函数)(x P ，用)(x P 来近似)(x f .这就是插值所要研究的问题.)(x P 称为)(x f 的插值函数.常用的插值函数是代数多项式或分段代数多项式.1.1 Lagrange 插值 1.1.1 方法介绍Lagrange 插值方法即，给定n 个插值节点以及对应的函数值信息，)(i i x f y =),2,1,0(n i =，利用n 次Lagrange 插值多项式公式，则对插值区间内任意x 的函数值y 可通过下式近似求得：)()(11∏∑≠==--=nkj j jk j n k k x x x x y x y .其中∏≠=--nkj j jkjx xx x 1称为插值基函数.可见，在Lagrange 插值中，对应1+n 个节点的插值基函数一共有1+n 个，每个基函数是一个n 次多项式. 1.1.2 MATLAB 实现 Lagrange.m% y=lagrange(x0,y0,x)，表示对x0，y0的节点信息使用Lagrange插值所得的多项式在x点的取值.function y=lagrange(x0,y0,x)n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end>> x=[0.4:0.1:0.8];>> clear>> x=[0.4:0.1:0.8];>> y=[-0.9 -0.6 -0.5 -0.35 -0.22];>> lagrange(x,y,0.54)ans =-0.55721.2分段低次插值1.2.1 方法介绍9596 利用给出的节点信息使用多项式插值，一般认为插值多项式的次数越高，对原函数的近似越好，但事实并非如此.19世纪Runge给出了一个反例，说明了插值多项式次数高并不意味着近似程度好.为了避免误差大的缺点，通常采用分段低次插值的做法.所谓分段低次插值，就是先把插值区间分成若干段，在每一段上进行低次多项式插值.1.2.2 MATLAB实现MATLAB自身提供了用于插值的内部函数interp1.该命令的使用方法为y=interp1(x,y,xi,’method’),其中(x,y)是已知的一组节点信息，xi表示插值函数在xi点取值，method表示所使用的插值方法，可选的方法有：● nearest 线性最近项插值；● linear 线性插值；● spline 三次样条插值；● cubic 三次插值.例：正弦函数的插值.>> x=[0:0.1:10];>> y=sin(x);>> xi=0:.25:10;>> yi=interp1(x,y,xi);>> plot(x,y,'o',xi,yi)971.3 Hermite 插值 1.3.1 方法介绍不少实际问题中不但要求在节点上函数值相等，还要求导数值也相等，甚至要求高阶导数值也相等，满足这一要求的插值多项式是Hermite 插值多项式.以下只给出函数值与一阶导数值个数相等且已知的情形.已知n 个插值节点n x x x ,,,21 及其对应的函数值n y y y ,,,21 和一阶导数值''2'1,,,ny y y .则计算插值区域内任意x 的函数值y 的Hermite 插值公式为 ∑=+--=ni i i i i i i y y y a x x h x y 1'])2)([()(,其中∑≠=≠=-=--=nji i ji i nij j ji j i x x a x x x x h 1211,)(.1.3.2 MATLAB 实现 hermite.m% 使用方法 y=hermite(x0,y0,y1,x),其中x0表示已知节点，% y0表示已知节点上的函数值信息，y1表示已知节点上的函数导数值信息, % x 表示插值多项式在x 上取值. function y=hermite(x0,y0,y1,x) n=length(x0),m=length(x); for k=1:m yy=0.0; for i=1:n h=1.0; a=0.0; for j=1:n if j~=i98 h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2; a=1/(x0(i)-x0(j))+a; end endyy=yy+h*((x(i)-x0(k))*(2*a*y0(i)-y1(i))+y0(i)); end y(k)=yy; end>> x0=[0.3 0.32 0.35]; >> y0=[0.3 0.3 0.34]; >> y1=[0.9 0.94 0.93]; >> x=[0.3:0.005:0.35]; >> y=hermite(x0,y0,y1,x);1.4 三次样条插值 1.4.1 方法介绍前面讨论的分段低次插值函数有一致收敛性，但光滑性较差.样条函数可以给出光滑的插值函数.其中重要的是三次样条函数.设区间],[b a 上给定有关划分,10b x x x a n =<<<= S 为],[b a 上满足下面条件的函数： ●),(2b a C S ∈；●S 在每个子区间],[1+i i x x 上是三次多项式. 则称S 为关于划分的三次样条函数.要确定三次样条函数，还需要边界条件.常用的边界条件有三种类型： ●);(')('),(')('00n n x f x S x f x S == ●);('')(''),('')(''00n n x f x S x f x S ==99●.,,1,0),()(0n j x S x S n j j ==1.4.2 MATLAB 实现interp1函数中，method 参数选为’spline ’即可，或者使用MATLAB 自带的spline 函数.下面举例说明spline 函数的应用. >>x=0:12; >>y=tan(pi*x/25); >>xi=linspace(0,12); >>yi=spline(x,y,xi);此例利用数据x 和y ，通过三次样条插值，计算得到正切函数的插值多项式在xi 处的取值. 1.2数值积分与数值微分许多实际问题需要计算函数的积分和微分，理论上计算积分的方法是利用微积分基本定理，计算微分的方法是利用微分的定义.但在实际应用中，函数的表达式往往是未知的，或者即使函数表达式已知，也未必满足理论公式的条件，并且在计算机实现上，还需要把连续的问题离散化.因此，我们要使用数值的方法来计算积分和微分.数值积分和数值微分的基本思想是，利用已有的节点以及节点上函数值的信息，先用插值多项式近似未知函数，再用我们熟知的方法求插值多项式的积分和微分，作为未知函数的积分和微分的近似.2.1 数值积分 2.1.1 方法介绍 对于积分⎰badx x f )((这里只考虑一元函数的积分)，将积分区间],[b a 划分为n等分，步长nab h -=，选取节点kh a x k +=，),,2,1,0(n k =，利用前面提到的多项式插值理论，可构造插值型的求积公式100 ∑=-=nk k n k n x f C a b I 0)()()(称为Newton －Cotes 公式，其中)(n k C 称为Cotes 系数，只与k n ,有关，可以通过计算得到.下表给出常用的几个Cotes 系数：特别地，Newton －Cotes 公式 ●当1=n 时，称为矩形公式； ●当2=n 时，称为梯形公式； ●当3=n 时，称为Simpson 公式； ●当4=n 时，称为Cotes 公式.如果直接应用上述公式，计算结果的误差还比较大，原因是插值的节点个数较少，在积分区间较大的情况下会产生很大的误差.实际应用中的做法是先将积分区间作一定的分划，然后在每个区间上应用上述公式，这样可以提高精度.改进后的公式分别称为复化矩形公式、复化梯形公式、复化Simpson 公式和复化Cotes 公式.另外,为了得到更高的代数精度,还有Gauss 求积公式和Romberg 求积公式,细节可参考相关文献资料.2.1.2 MATLAB 实现MATLAB 中自带了求积分的几个函数,这里我们把它们列表如下：101例：求积分⎰+124dx x x. fun.mfunction y=fun(x) y=x./(4+x.^2); >>quad(’fun ’,0,1) ans = 0.1116除数值积分函数外，MATLAB 还提供了两个符号积分函数int 和symsum 。
关于它们的使用方法可参看MATLAB 的帮助文件，这里不做展开.2.2 数值微分微分和差分也是进行数学计算的重要部分，但由于计算中很少出现稳定性和精度问题，所以相应的算法也较少。
MATLAB 提供了几个功能函数用于解决微分和和差分问题。
下面举例说明diff 和jacobian 这两个函数的用法： ● diff例：求)sin(2x 的导数. >> x=sym('x'); >> diff(sin(x^2)) ans = 2*cos(x^2)*x ● jacobian例：求如下函数的Jacobi 矩阵.⎪⎪⎪⎭⎫⎝⎛+-++=z e z x e y x z y x F y z cos 5sin 2),,(2>> x=sym('x');y=sym('y');z=sym('z');>> f1='[2*x+sin(y)+exp(z);x^2-5*z;exp(y)+cos(z)]';102 >> jacobian(f1,[x,y,z]) ans =[ 2, cos(y), exp(z)] [ 2*x, 0, -5] [ 0, exp(y), -sin(z)] 1.3求解非线性方程组非线性方程的求根一般采用遍历法、二分法和迭代法求解.遍历法由于计算量大，精度差，一般只有搜索到根区间时使用.下面主要介绍二分法和迭代法，迭代法主要介绍不动点迭代和Newton 迭代.3.1 二分法二分法的理论依据是闭区间上连续函数的介值定理，即：如果)(x f 是],[b a 上的连续函数，且0)()(<b f a f ，那么在],[b a 至少有一个根.对此区间二分，再用介值定理找到有根区间，再二分，直到有根区间长度足够小，二分结束，用当前区间的中点作为根的近似.这就是二分法.由于该算法较为简单，在计算机上的实现这里不做介绍. 3.2 不动点迭代 3.2.1 方法介绍设一元函数)(x f 是连续的，要求解0)(=x f ，为了进行迭代，变换方程形式为)(x x ϕ=，于是构造迭代公式)(1k k x x ϕ=+.如果)(*∞→→k x x k ，则称此迭代为不动点迭代.但此迭代是否收敛与)(x ϕ的形式有关，并且如果收敛，收敛速度也受)(x ϕ影响.3.2.2 MATLAB 实现 g.mfunction y=g(x) y=log(3*x.^2);103iterate.m% y=iterate(x)，其中x 是迭代的初值 function y=iterate(x) x1=g(x); n=1;while (abs(x1-x)>=1.0e-6)&(n<=1000) x=x1; x1=g(x); n=n+1; end y=x1; >>iterate(3) ans= 3.73313.2 Newton 迭代 3.2.1 方法介绍Newton 法是最重要，而且是应用最为广泛的一种迭代法.在这里对其推导不做介绍.Newton 法的迭代公式为：)(')(1k k k k x f x f x x -=+.Newton 法是局部平方收敛的，因此使用它计算更快捷. 3.2.2 MATLAB 实现类似不动点迭代，可以编写Newton 迭代的MATLAB 程序，这里还要计算的一项是导函数.代码的细节略去.另外，在MATLAB 的Symbolic Toolbox 中提供了用于求解非线性方程（组）的函数fsolve ，使用极其简便.例：fc.mfunction y=fc(x)y(1)=x(1)-0.7*sin(x(1))-0.2cos(x(2));y(2)=x(2)-0.7*cos(x(1))+0.2sin(x(2));y=[y(1) y(2)];>> x0=[0.5 0.5];>> fsolve('fc',x0)ans =0.5265 0.5079§2数值代数求解线性方程组和求解矩阵的特征值问题的算法很多，细节可以参考数值分析的书籍.而在MATLAB中，有着现成的函数可供调用，这些函数集许多优秀算法于一身，而且是算法自适应的，即根据不同的已知条件选择合适的算法来计算.2.1求解线性方程组Ax 的解法一般可分为两类：一是直接法，通过矩阵的变形、关于线性方程组b消去直接得到方程的解，这类方法是解低阶稠密矩阵方程组的有效方法；二是迭代法，就是用某种极限过程去逐渐逼近方程组精确解的方法，迭代法是解大型稀疏矩阵方程组的重要方法.1.1 直接法关于线性方程组的直接解法，常见的有Gauss消去法、选主元消去法、平方根法.追赶法等等，大多数数值分析的书都有比较详尽的介绍.在MATLAB中，只须用一个“\”就解决问题.虽然表面上是一个简单的符号，但内部却包含了许多算法，如对超定方程组使用最小二乘法，对不定方程组它将给出范数最小的一个解，解三对角方程组使用追赶法等等.具体用法为：初始化矩阵A和右端向量b，那么>>x=A\b1.2 迭代法104105常见的迭代法有Jacobi 迭代、Gauss －Seidel 迭代和SOR 迭代，这里我们只介绍Jacobi 迭代法，另外两个迭代法与Jacobi 迭代法类似，可在相关参考书上找到算法的细节.设A 是n 阶非奇异方阵，把A 写成U L D A --=，其中D 是对角矩阵，它的对角元对应A 的对角元，而L 和U 分别是A 的严格下三角部分和严格上三角部分（不包括对角元），则b D x U L D x 11)(--++=，由此构造迭代法f Bx x k k +=+1，其中)(1U L D B +=-，b D f 1-=.Jacobi 迭代法的MATLAB 实现如下：jacobi.m% 使用方法y=jacobi(a,b,x0)% 其中a 是系数矩阵，b 是右端向量，x0是迭代的初始向量function y=jacobi(a,b,x0)D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;while norm(y-x0))>=1.0e-6x0=y;y=B*x0+f;n=n+1;endy106 例：求解方程组⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛----679102021010110321x x x>>a=[10 -1 0;-1 10 -2;0 -2 10];>>b=[9;7;6];>>jaboci(a,b,[0;0;0])y=0.99580.95790.79162.2 特征值问题在实际应用中，求解特征值问题的常见算法有幂法、QR 算法等，其中幂法及其变形（反幂法）是求解矩阵端部特征值和对应特征向量的有效方法，QR 算法是求解矩阵所有特征值和特征向量的迭代法.在MATLAB 中，主要应用eig 函数或eigs 函数（针对稀疏矩阵）来求解特征值问题.例：>>A=rand(3);>>[V,d]=eig(A)V =-0.6571 -0.7865 0.7614-0.2275 0.3771 -0.6380-0.7186 0.4892 0.1154d =1.6175 0 00 0.4332 00 0 0.6121107 上例计算随机生成的三阶矩阵的特征值和特征向量，返回V 表示特征向量构成的矩阵，D 表示以所有特征值为对角元的对角矩阵.如果求解线性方程组或特征值问题中的矩阵是大型稀疏的，有专门的算法可供选择，例如投影类的子空间迭代法，算法的实现可调用现有的一些软件包.§3 微分方程数值解科学技术中常常需要求解常微分方程的定解问题，这类问题最简单的形式是下面的一阶方程的初值问题：⎩⎨⎧==.)(),,('00y x y y x f y 这一部分我们介绍其数值解法以及在MATLAB 中的应用.所谓数值解法，就是寻求)(x y 在一系列离散节点<<<<<<+1321n n x x x x x上的近似值 ,,,,,121+n n y y y y .这里我们介绍Euler 方法和Runge －Kutta 方法.1、 Euler 方法对初值问题进行变换，用差商近似导数，有),(11n n nn n n y x f x x y y =--++， 进而得到),(1n n n n y x hf y y +=+，这里h 表示相邻两个节点1,+n n x x 之间的距离.这就是著名的Euler 格式.如果给定初值0y ，那么可以逐步算出 ,,21y y .Euler 格式实际上是做了线性近似.2、 Runge －Kutta 方法2.1 方法介绍Runge －Kutta 方法实质上是间接地使用Taylor 展开二阶近似的一种方法（Taylor 级数展开到二次项），如果解的光滑性较好，那Runge －Kutta 方法的精度比Euler 法要高.108 经典的四阶Runge－Kutta格式为⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧++=++=++==++++=+).,(),2,2(),2,2(),,(),22(6342312142211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn但需要注意的是，Runge－Kutta方法基于函数的Taylor展开，因此要求所求的解具有很好的光滑性，反之，如果解的光滑性差，那么使用Runge－Kutta方法可能效果会比较差.在实际计算中，应针对问题的具体特点选择合适的算法.2.2 MATLAB实现MATLAB中有几个专门用于解常微分方程的功能函数，如ode23、ode45、ode23s 等.ode23系列采用二阶、三阶的Runge－Kutta方法，ode45系列则采用四阶、五阶的Runge－Kutta方法.以ode45为例演示该函数的用法：fun.mfunction f=fun(x,y)f=-2*y+2*x.^2+2*x;>> ode45('fun',[0,0.5],1)。
