第五章数值计算.doc

《第五章数值计算.doc》由会员分享，可在线阅读，更多相关《第五章数值计算.doc（17页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第五章 数值计算5.1多项式的计算在全部初等函数里，多项式是结构最为简单，应用最为广泛的一类函数，有必要单独列出进行讨论。1多项式的表示法在MATLAB中将多项式采用降幂排列，并可表示为向量的形式。例如，多项式p(x)=x3-12x+5可用向量表示为p(x)=1 0 12 52多项式在给定点的值利用MATLAB中的polyval函数，可以直接求出多项式在给定点的值。y=polyval(p,x0),其中p为多项式，x0为给定点。【例51】求多项式p=4 -3 2 1在点x=1 2 3 4 5处的值输入p=4 -3 2 1;x=1 2 3 4 5y=polyval(p,x); %求多项式在x处的值

2、。x,y运行后可得x = 1 2 3 4 5y = 4 25 88 217 4363多项式方程求根利用MATLAB中的roots函数，可以直接求得多项式方程的根。p(xi)=0的根xi用函数xi=roots(p).【例52】求多项式方程的根输入p1=1 0 -12 5;p2=1 0 -12 20;x1=roots(p1);x2=roots(p2);x1,x2运行后可得x1 = -3.6562 3.2332 0.4230x2 = -4.1072 2.0536 + 0.8075i 2.0536 - 0.8075i4多项式相乘两多项式相乘，可以利用MATLAB中的conv(convolution)函

3、数来求得.p3=conv(p1,p2),5、多项式相除：两多项式与相除，可以利用MATLAB中的deconv 函数来求得，其调用格式为q,r=deconv(p1,p2),其中p1为被除式，p2为除式，q为商式，r为余式。【例53】多项式相乘与相除输入p1=5 -4 3 -2 1 -6;p2=1 0 -2 3;p3=conv(p1,p2);q,r=deconv(p1,p2);p3,q,r运行后可得p3 = 5 -4 -7 21 -17 7 -8 15 -18q = 5 -4 13r = 0 0 0 -25 39 -456有理分式的分解与合并：有理分式同多项式有着紧密的联系，依照代数定理，任何有理

4、分式都可分解为部分分式（或称最简分式）之和式中pm(x)与qn(x)为m次与n次多项式，当mn时，r(x)为商式，当mn时，r(x)为空，bi为方程qn(x)=0的根。ai为对应的实数或复数。将有理分式进行分解，是一件十分复杂和繁琐的工作。如果我们调用MATLAB里的residule函数，将会使问题变得简单容易，其调用格式为a,b,r=residule(p,q); 有理分式分解.p,q=residule(a,b,r). 将简单分式之和合并为有理分式.【例54】将有理分式分解为最简分式之和。输入p=1 0 0 0;q=1 1 -2; a,b,r=residue(p,q);a,b,r运行后可得a

5、= 2.6667 0.3333b = -2 1r = 1 -1即：【例55】将有理分式分解为最简分式之和。输入p=1 1;q=1 -6 11 -6;a,b,r=residue(p,q);a,b,r运行后可得a = 2.0000 -3.0000 1.0000b = 3.0000 2.0000 1.0000r = 故求得【例56】将简单分式之和合并为有理分式。输入a=2.667 0.333;b=-2 1;r=1 -1;p,q=residue(a,b,r);p,q运行后可得p = 1.0000 0 0 -0.0010q = 1 1 -252导数的数值计算导数的数值计算或称数值求导，是指导数值的一类近

6、似计算方法，它不必去求出导函数（需要符号求导运算），而是利用函数值的差分来构造出近似导数值。下面为求导的几个常见命令：1、yx=diff(y)./diff(x)（即y=y/x）其中diff=differential，差别的，微分。2、fx,fy=gradient(f, x, y)Gradient,梯度，斜坡。f为已给的二元函数， x, y 为自变量x与y的步长，fx,fy为二元函数的两个偏导数的近似值。3、quiver(fx,fy)是二维向量场的表现函数。qurver,大群，箭袋。【例57】导数的数值计算。输入x=0:0.1:3y=x.2yx=gradient(y,0.01)yx0=diff(

7、y)./diff(x)x=0:0.5:3;y=x.2;yx=diff(y)./diff(x);x,y,yx运行后可得x = 0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000y = 0 0.2500 1.0000 2.2500 4.0000 6.2500 9.0000yx = 0.5000 1.5000 2.5000 3.5000 4.5000 5.5000【例58】画函数f=x.2+x.*y的二维向量场输入x,y=meshgrid(0:0.2:1,0:0.3:1.5);f=x.2+x.*y; fx,fy=gradient(f,0.2,0.3);fx,fyc

8、ontour(f)hold on %保留当前图形。quiver(fx,fy)hold off运行后结果如图5-1图5-1【例59】画函数z=x.*exp(-x.2-y.2)的二维向量场输入 x,y=meshgrid(-2:.2:2); z=x.*exp(-x.2-y.2); zx,zy=gradient(z,0.2,0.2);figure(1)mesh(x,y,z)figure(2)contour(z);hold on quiver(zx,zy) zx,zy 图5-2如图5-25.3定积分的数值计算对于定积分，如果相应的不定积分不易求出，或者根本不能表示为初等函数时，那么就只能用近似数值计算的

9、办法去求定积分的值了。函数 quad、quadl、quad8功能 数值定积分，自适应Simpleson积分法。格式 q = quad(fun,a,b,tol) %用指定的绝对误差tol代替缺省误差。Tol越大，函数计算的次数越少，速度越快，但结果精度变小。q,n = quad(fun,a,b,) %同时返回函数计算的次数n = quadl(fun,a,b,) %用高精度进行计算，效率可能比quad更好。【例510】求y=sin(x)./x的积分。输入function y=f(x)y=sin(x)./x;s=quad(p99l1,1,10)%应在工作窗口求解。或在另一个editor/debugg

10、er中求解。运行后可得s = 0.7123【例511】求y= exp(-x.2/2)的积分。输入f=inline(exp(-x.2/2),x);s=quad(f,-8,8)运行后可得s = 2.5066【例512】求y=3*x.2./(x.3-2*x.2+3的积分。输入fun = inline(3*x.2./(x.3-2*x.2+3);Q1 = quad(fun,0,2)Q2 = quadl(fun,0,2)计算结果为：Q1 = 3.7224Q2 = 3.7224同定积分的概念一样，二重积分的值不只依赖于被积函数，同时也要依赖于积分区域的情况，如果为矩形区域，则计算过程将会变得十分简单，如果为

11、一般区域，则计算过程将会比较复杂，这里不进行讨论。函数 dblquad功能 矩形区域上的二重积分的数值计算格式 q = dblquad(fun,xmin,xmax,ymin,ymax) %调用函数quad在区域xmin,xmax, ymin,ymax上计算二元函数z=f(x,y)的二重积分。输入向量x，标量y，则f(x,y)必须返回一用于积分的向量。【例513】求 y./sin(x)+x.*exp(y)的二重积分。输入fun = inline(y./sin(x)+x.*exp(y);Q = dblquad(fun,1,3,5,7)计算结果为：Q = 3.8319e+0035.4非线性方程的数值

12、解如果非线性方程（组）是多项式形式，求这样方程（组）的数值解，可以直接调用本章5.2节中已介绍过的函数，如果非线性方程（组）中出现有超越函数，例如等，则无法使用函数，而必须采用系统中提供的另一函数来求解。函数当然也可用于多项式方程（组）的求解，但它的计算量要比函数明显地多。【例514】求方程x3-12*x+5=0的一个数值解。建立M文件如下：function f=fun(x) %fun为函数名，也是保存的文件名f=x3-12*x+5在命令窗口输入：x0=1x1=fsolve(fun,x0)f=fun(x1) %根满足方程的情况计算结果为：x1 = 0.4230f = 4.4409e-015【例

13、515】求方程组x3-y2=0;exp(-x)-y=0的一组数值解。建立M文件如下：function f=fun(t) %fun为函数名，也是保存的文件名x=t(1);y=t(2);f(1)=x3-y2;f(2)=exp(-x)-y;在命令窗口输入：t0=1,1t=fsolve(fun,t0)f=fun(t) %根满足方程的情况计算结果为：t = 0.6488 0.5226f = 1.0e-006 * 0.2196 0.0321【例516】求方程组x+y+z-6=0;x+y*z+z*x-8=0;exp(-x)+log(y)+z-2=0;的一组数值解。建立M文件如下：function f=fun

14、(t) %fun为函数名，也是保存的文件名x=t(1);y=t(2);z=t(3);f(1)=x+y+z-6;f(2)=x+y*z+z*x-8;f(3)=exp(-x)+log(y)+z-2;在命令窗口输入：t0=1,1 ,1t=fsolve(fun,t0)f=fun(t) %根满足方程的情况计算结果为：t = 2.6132 2.2877 1.0992f = 1.0e-012 * 0.0036 -0.1883 -0.08665.5级数的数值和对于有限项和与无穷项和，当通项中含有任意参数或时，此种含有符号变量的求法我们已在前面进行了讨论。当通项中不含有任意参数或时，此种和式的结果将是一个数值，可

15、称之为级数的数值和。【例517】求级数fn=1/(n*n)的数值和。syms n;fn=1/(n*n);s10=symsum(fn,n,1,10);sinf=symsum(fn,n,1,inf);s10,sinfs10 =/sinf =1/6*pi25.6 一元插值在生产实践和科技活动中，人们常常从实验或测量中行到数据，称为原始数据，而将原始数据所确定的点列，自负盈亏为原始点列，或型值点列。人们希望将这些数据或点列中所蕴含的客观规律用一个函数或者一条平面曲线表达出来，这样的函数或曲线称为理想函数或理想曲线。插值法是实用的数值方法，是函数逼近的重要方法。在生产和科学实验中，自变量x与因变量y的函

16、数y = f(x)的关系式有时不能直接写出表达式，而只能得到函数在若干个点的函数值或导数值。当要求知道观测点之外的函数值时，需要估计函数值在该点的值。如何根据观测点的值，构造一个比较简单的函数y=(x)，使函数在观测点的值等于已知的数值或导数值。用简单函数y=(x)在点x处的值来估计未知函数y=f(x)在x点的值。寻找这样的函数(x)，办法是很多的。(x)可以是一个代数多项式，或是三角多项式，也可以是有理分式；(x)可以是任意光滑（任意阶导数连续）的函数或是分段函数。函数类的不同，自然地有不同的逼近效果。在许多应用中，通常要用一个解析函数（一、二元函数）来描述观测数据。根据测量数据的类型：1测

17、量值是准确的，没有误差。2测量值与真实值有误差。这时对应地有两种处理观测数据方法：假定数据测量是精确时，一般用插值法，否则用曲线拟合。一元多项式插值的概念是：如何在整段区间上构造一个多项式函数；或者在分段区间上，构造次数较低的多项式函数。已知原始数据（离散数据）x=x1,x2,x n+1，y=y1,y2y n+1所确定的型值点列M1(x1,y1),M2(x2,y2)M n+1(x n+1,y n+1)，我们要用它构造一条“最”光滑的曲线，如果曲线通过所有的型值点，称为插值；如果要求曲线逼近型值点称为拟合。（因为 型值点本身就可能带有各种误差，所以要求曲线通过所有的型值点列有时是无必要的。）插值

18、与拟合的区别在于：1、插值通过所有的型值点；拟合逼近型值点。2、插值与拟合都是用多项式构造光滑的理想曲线，在插值时，不同型值点列间的多项式不同，而在拟合时，整个区间上都用一个多项式，因些，拟合曲线可用多项式表示。插值命令命令 interp1功能 一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。格式 yi = interp1(x,Y,xi,method) %用指定的算法计算插值：nearest：最近邻点插值，直接完成计算；linear：线性插值（缺省方式），直接完成计算；spline：三次样条函数插值。命令生成一系列用

19、于分段多项式操作的函数。命令spline用它们执行三次样条函数插值；【例518】一元插值例。x0=0:0.2:1;y0=0 0.32 0.21 0.1 0.04 0.02;x=0:0.02:1;y1=interp1(x0,y0,x);y2=interp1(x0,y0,x,cubic);y3=interp1(x0,y0,x,spline);y4=interp1(x0,y0,x,nearest);plot(x0,y0,ko,x,y1,b:,x,y2,r-,x,y3,m-,x,y4,b-)legend(origin data,linear,cubic,spline,nearest);运行后结果如图5

20、-3 图5-3【例519】一元插值例。x0=0:1:10; y0=cos(x0);x=0:0.1:10;y1=interp1(x0,y0,x);y2=interp1(x0,y0,x,cubic);y3=interp1(x0,y0,x,spline);y4=interp1(x0,y0,x,nearest); 图5-4plot(x0,y0,O,x,y1,k-,x,y2,b:,x,y3,r-,x,y4,c);legend(origin data,linear,cubic,spline,nearest);如图5-4【例520】一元插值例。x0=0:0.4:4;y0=0 0.1 0.85 1.05 3.

21、08 3.5 6.07 7.9 9.7 13.96 14.68;x=0:0.04:4;y=interp1(x0,y0,x,spline);plot(x0,y0,o,x,y);legend(origin data,spline)运行后结果如图5-5图5-55.7一元拟合多项式插值要求多项式曲线必须穿过所讨论区间上每个型值点，如果在全区间上插值，当多项式的次数较大时常会引起曲线在区间端点与附近的剧烈振荡（称为runge现象），如果在子区间上插值，每个子区间上多项式的次数虽然降低了，但在连接点处曲线的光滑性也随之降低，这种矛盾常常无法得到统一的解决。从实验者的角度看，要求逼近曲线穿过每个型值点，既无

22、此种必要，又无进一步提高逼近精度的可能，因为型值点本身就带有各种误差。人们想到，能否构造一条次数较低的多项式曲线，只要求尽可能地贴近全部型值点，而不必每点穿过。“尽可能地贴近”通常选用最小二乘意义下的逼近，这就是一元拟合的概念。MATLAB为用户提供了一元拟合具体实现的调用函数，它的调用格式如下：p=polyfit(x0,y0,r)p,s=polyfit(x0,y0,r)式中，向量x0与y0的分量均为已知，即前面所说的型值点列；r为拟合多项式的估计次数；s为估计误差用的参数；p为输出向量，其分量为求得的r次多项式的系数（按降幂排列）【例521】拟合原理。x0=-4,-2,0,2,4 y0=3.

23、*x0.2+5.*x0y2=30.5,0,3,19.5,73x=-5:0.1:5;y1=interp1(x0,y0,x,spline);p1=polyfit(x0,y0,1);y11=polyval(p1,x);y3=polyval(p1,x0)plot(x0,y2,mo,x0,y3,p,x,y11,x0,y0,h,x,y1) xlabel(x)ylabel(y)legend(origin data,x0y3,x y,x0 y0,x y1) 图5-6如图5-6【例522】一元拟合例x0=0:0.4:4;y0=0 0.1 0.85 1.05 3.08 3.5 6.07 7.9 9.7 13.96

24、 14.68;p1=polyfit(x0,y0,1);p2=polyfit(x0,y0,2); p3=polyfit(x0,y0,3);p10=polyfit(x0,y0,10);p1,p2,p3,p10x=0:0.1:4;y1=polyval(p1,x);y2=polyval(p2,x);y3=polyval(p3,x);y10=polyval(p10,x);plot(x0,y0,ko,x,y1,b:,x,y2,r-,x,y3,k-,x,y10,k-);legend(y0-x0,y1-x,y2-x,y3-x,y10-x)5.8二元函数的插值一元函数插值的概念可以完全类似地推广到二元函数上去，

25、只是型值点给出的形式更加复杂一些。二元插值的任务是：在矩形域上，或者在的分片子域上构造二元多项式去逼近由点列所蕴含的理想函数，或说构造多项式曲面去逼近理想曲面。二元多项式通常取为双三次的。它的使用格式如下：Zij=interp2(x0,y0,z0,xi,yj,插值方法)2表示二元函数的插值。X0,y0,z0为已知的型值点。Xi,yj为插值点。Zij为二元多项式在插值点处的函数值。插值方法有：linear：双线性插值算法（缺省算法）；nearest：最临近插值；spline：三次样条插值；cubic：双三次插值。【例523】二元插值例。x0=-3:1:3;y0=-3.6:1.2:3.6;z0=0

26、.55 1.25 0.75 -0.35 -1.45 -1.95 -1.25; 2.4 3.1 2.6 1.5 0.4 -0.1 0.6; 2.17 2.87 2.37 1.27 0.17 -0.33 0.37; 0.9 1.6 1.1 0 -1.1 -1.6 -0.9; -0.37 0.33 -0.17 -1.27 -2.37 -2.87 -2.17; -0.6 0.1 -0.4 -1.5 -2.6 -3.1 -2.4; 1.25 1.95 1.45 0.35 -0.75 -1.25 -0.55;x,y=meshgrid(-3:0.3:3,-3.6:0.36:3.6);z1=interp2(x

27、0,y0,z0,x,y);z2=interp2(x0,y0,z0,x,y,cubic);z3=interp2(x0,y0,z0,x,y,spline);figure(1);plot3(x0,y0,z0,mo);figure(2);mesh(x,y,z1);figure(3);mesh(x,y,z2);figure(4);mesh(x,y,z3);figure(5);subplot(2,2,1)plot3(x0,y0,z0,mo);legend( origin data )grid onsubplot(2,2,2)mesh(x,y,z1); 图5-7legend(linear)subplot(2

28、,2,3)mesh(x,y,z2);legend(cubic)subplot(2,2,4)mesh(x,y,z3);legend(spline) 运行后结果如图5-7，图5-8 图5-8【例524】二元插值例。x,y = meshgrid(-3:.25:3);z = peaks(x,y);xi,yi = meshgrid(-3:.125:3);ZZ = interp2(x,y,z,xi,yi);mesh(x,y,z);hold on;mesh(xi,yi,ZZ+15)axis(-3 3 -3 3 -5 20);shading flathold off运行后结果如图5-9图5-9【例525】二元插值例。years = 1950:10:1990;service = 10:10:30;wage = 150.697 199.592 187.625 179.323 195.072 250.287 203.212 179.092 322.767 226.505 153.706 426.730 249.633 120.281 598.243;w = interp2(service,years,wage,15,1975)w = 190.6288 （白洪波 孙福玉）