《数值计算方法精讲.doc》由会员分享，可在线阅读，更多相关《数值计算方法精讲.doc（35页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、【精品文档】如有侵权，请联系网站删除，仅供学习与交流数值计算方法精讲.精品文档.数值计算方法5.1引言本章将花较大的篇幅讨论若干常见数值计算问题：线性分析、一元和多元函数分析、微积分、数据分析、以及常微分方程求解等。但与一般数值计算教科书不同，本章的讨论重点是：如何利用现有的世界顶级数值计算资源MATLAB。至于数学描述，本章将遵循“最低限度自封闭”的原则处理，以最简明的方式阐述理论数学、数值数学和MATLAB计算指令之间的内在联系及区别。对于那些熟悉其他高级语言（如FORTRAN，Pascal，C+）的读者来说，通过本章，MATLAB卓越的数组处理能力、浩瀚而灵活的M函数指令、丰富而友善的图
2、形显示指令将使他们体验到解题视野的豁然开朗，感受到摆脱烦琐编程后的眉眼舒展。对于那些经过大学基本数学教程的读者来说，通过本章，MATLAB精良完善的计算指令，自然易读的程序将使他们感悟“教程”数学的基础地位和局限性，看到从“理想化”简单算例通向科学研究和工程设计实际问题的一条途径。对于那些熟悉MATLAB基本指令的读者来说，通过本章，围绕基本数值问题展开的内容将使他们体会到各别指令的运用场合和内在关系，获得综合运用不同指令解决具体问题的思路和借鉴。由于MATLAB的基本运算单元是数组，所以本章内容将从矩阵分析、线性代数的数值计算开始。然后再介绍函数零点、极值的求取，数值微积分，数理统计和分析，
3、拟合和插值，Fourier分析，和一般常微分方程初值问题。本章的最后讨论稀疏矩阵的处理，因为这只有在大型问题中，才须特别处理。从总体上讲，本章各节之间没有依从关系，即读者没有必要从头到尾系统阅读本章内容。读者完全可以根据需要阅读有关节次。除特别说明外，每节中的例题指令是独立完整的，因此读者可以很容易地在自己机器上实践。5.1 LU分解和恰定方程组的解5.1.1 LU分解、行列式和逆（1）LU分解 （2）行列式和逆5.1.2 恰定方程组的解【*例5.2.2-1】“求逆”法和“左除”法解恰定方程的性能对比（1）为对比这两种方法的性能，先用以下指令构造一个条件数很大的高阶恰定方程。rand(stat
4、e,12);%选定随机种子，目的是可重复产生随机阵A。A=rand(100,100)+1.e8;%rand(100,100)生成（100100）均匀分布随机矩阵。%每个随机阵元素加的目的是使A阵条件数升高。x=ones(100,1);%令解向量 x 为全1的100元列向量。b=A*x;%为使 Ax=b 方程一致，用A和 x 生成 b 向量。cond(A) %求A阵的条件数。 ans = 1.4426e+012 （2）“求逆”法解恰定方程的误差、残差、运算次数和所用时间flops(0);tic%浮点运算计数器置0 ；启动计时器Stopwatch Timerxi=inv(A)*b;% xi 是用“
5、求逆”法解恰定方程所得的解。ti=toc%关闭计时器，并显示解方程所用的时间。ci=flops%“求逆”法解方程所用的运算次数eri=norm(x-xi)%解向量 xi 与真解向量 x 的范-2误差。rei=norm(A*xi-b)/norm(b)%方程的范-2相对残差 ti = 0.9300ci = 2070322eri = 3.0708e-004rei = 6.6280e-007 （3）“左除”法解恰定方程的误差、残差、运算次数和所用时间flops(0);tic;xd=Ab;%是用“左除”法解恰定方程所得的解。td=toc,cd=flops,erd=norm(x-xd),red=norm(
6、A*xd-b)/norm(b) td = 0.2200cd = 741872erd = 3.2243e-004red = 2.0095e-016 5.1.3 范数、条件数和方程解的精度【*例5.2.3-1】Hilbert矩阵是著名的病态矩阵。MATLAB中有专门的Hilbert矩阵及其准确逆矩阵的生成函数。本例将对方程近似解和准确解进行比较。所谓n阶Hilbert矩阵的形式是： 。N=6 8 10 12 14;%本例计算的矩阵阶数for k=1:length(N)n=N(k);%矩阵的阶H=hilb(n);%产生n阶Hilbert矩阵Hi=invhilb(n);%产生完全准确的n阶逆Hilbe
7、rt矩阵b=ones(n,1);%生成n阶全1向量x_approx=Hb; %利用左除H求近似解x_exact=Hi*b;%利用准确逆Hilbert矩阵求准确解ndb=norm(H*x_approx-b);nb=norm(b);ndx=norm(x_approx - x_exact);nx=norm(x_approx);er_actual(k)=ndx/nx;%实际相对误差K=cond(H);%计算Hilbert矩阵的条件数er_approx(k)=K*eps;%最大可能的近似相对误差er_max(k)=K*ndb/nb; %最大可能的相对误差end disp(Hilbert矩阵阶数),dis
8、p(N)format short edisp(实际误差 er_actual),disp(er_actual),disp()disp(近似的最大可能误差 er_approx),disp(er_approx),disp()disp(最大可能误差 er_max),disp(er_max),disp() Hilbert矩阵阶数 6 8 10 12 14实际误差 er_actual 5.0339e-011 8.5981e-008 2.2819e-004 1.3381e-001 3.9641e+000近似的最大可能误差 er_approx 3.3198e-009 3.3879e-006 3.5583e-0
9、03 3.9259e+000 3.4573e+002最大可能误差 er_max 6.0095e-007 2.4531e-002 1.4094e+003 2.9206e+007 2.4178e+010 5.2 矩阵特征值和矩阵函数5.2.1 特征值和特征向量的求取【例5.3.1-1】简单实阵的特征值问题。A=1,-3;2,2/3;V,D=eig(A) V = -0.7728 + 0.0527i -0.7728 - 0.0527i 0 + 0.6325i 0 - 0.6325iD = 0.8333 + 2.4438i 0 0 0.8333 - 2.4438i 【*例5.3.1-2】本例演示：如矩阵
10、中有元素与截断误差相当时的特性值问题。A=3 -2 -0.9 2*eps -2 4 -1 -eps -eps/4 eps/2 -1 0 -0.5 -0.5 0.1 1 ;V1,D1=eig(A);ER1=A*V1-V1*D1V2,D2=eig(A,nobalance);ER2=A*V2-V2*D2 ER1 = 0 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 1.1227ER2 = 1.0e-015 * 0.4441 -0.2220 0
11、.1471 -0.2220 0 0.0555 -0.3629 0.2776 -0.0172 -0.0015 0.0066 0 0 -0.2220 -0.1110 0.1388 【*例5.3.1-3】指令eig与eigs的比较。rand(state,1),A=rand(100,100)-0.5;t0=clock;V,D=eig(A);T_full=etime(clock,t0)%指令eig的运作时间。options.tol=1e-8;%为eigs设定计算精度。options.disp=0;%使中间迭代结果不显示。t0=clock;v,d=eigs(A,1,lr,options);%计算最大实部特
12、征值和特征向量。T_part=etime(clock,t0)%指令eigs的运作时间。Dmr,k=max(real(diag(D);%在eig求得的全部特征值中找最大实部的那个。d,D(1,1) T_full = 2.8000T_part = 19.5500d = 3.0140 - 0.2555ians = 3.0140 + 0.2555i vk1=V(:,k+1);%与d相同的特征向量应是V的第k+1列。vk1=vk1/norm(vk1);v=v/norm(v);%向量长度归一。V_err=acos(norm(vk1*v)*180/pi%求复数向量之间的夹角（度）。D_err=abs(D(k
13、+1,k+1)-d)/abs(d)%求两个特征值间的相对误差。 V_err = 8.5377e-007D_err =1. 7098e-010 5.2.2 特征值问题的条件数【例5.3.2-1】矩阵的代数方程条件数和特征值条件数。B=eye(4,4);B(3,4)=1;Bformat short e,c_equ=cond(B),c_eig=condeig(B) B = 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 1c_equ = 2.6180e+000Warning: Matrix is close to singular or badly scaled. Results may
14、 be inaccurate. RCOND = 1.500000e-018. In E:MAT53toolboxmatlabmatfuncondeig.m at line 30c_eig = 1.0000e+000 1.0000e+000 3.3333e+017 3.3333e+017 【*例5.3.2-2】对亏损矩阵进行Jordan分解。A=gallery(5)%MATLAB设置的特殊矩阵，它具有五重特征值。VJ,DJ=jordan(A);%求出准确的特征值，使A*VJ=VJ*D成立。V,D,c_eig=condeig(A);c_equ=cond(A);DJ,D,c_eig,c_equ A
15、= -9 11 -21 63 -252 70 -69 141 -421 1684 -575 575 -1149 3451 -13801 3891 -3891 7782 -23345 93365 1024 -1024 2048 -6144 24572DJ = 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0D = Columns 1 through 4 -0.0328 + 0.0243i 0 0 0 0 -0.0328 - 0.0243i 0 0 0 0 0.0130 + 0.0379i 0 0 0 0 0.0130 - 0.0379i 0 0
16、0 0 Column 5 0 0 0 0 0.0396 c_eig = 1.0e+010 * 2.1016 2.1016 2.0251 2.0251 1.9796c_equ = 5.2133e+017 5.2.3 复数特征值对角阵与实数块特征值对角阵的转化【*例5.3.3-1】把例5.3.1-1中的复数特征值对角阵D转换成实数块对角阵，使VR*DR/VR=A。VR,DR=cdf2rdf(V,D) VR = -0.7728 0.0527 0 0.6325DR = 0.8333 2.4438 -2.4438 0.8333 5.2.4 矩阵的谱分解和矩阵函数【*例5.3.4-1】数组乘方与矩阵乘方的
17、比较。clear,A=1 2 3;4 5 6;7 8 9;A_Ap=A.0.3%数组乘方A_Mp=A0.3 %矩阵乘方 A_Ap = 1.0000 1.2311 1.3904 1.5157 1.6207 1.7118 1.7928 1.8661 1.9332A_Mp = 0.6962 + 0.6032i 0.4358 + 0.1636i 0.1755 - 0.2759i 0.6325 + 0.0666i 0.7309 + 0.0181i 0.8292 - 0.0305i 0.5688 - 0.4700i 1.0259 - 0.1275i 1.4830 + 0.2150i 【*例5.3.4- 2
18、】标量的数组乘方和矩阵乘方的比较。（A取自例5.3.4-1）pA_A=(0.3).A%标量的数组乘方pA_M=(0.3)A%标量的矩阵乘方 pA_A = 0.3000 0.0900 0.0270 0.0081 0.0024 0.0007 0.0002 0.0001 0.0000pA_M = 2.9342 0.4175 -1.0993 -0.0278 0.7495 -0.4731 -1.9898 -0.9184 1.1531 【*例5.3.4-3】sin的数组运算和矩阵运算比较。（A取自例5.3.4-1）A_sinA=sin(A)%数组运算A_sinM=funm(A,sin)%矩阵运算 A_si
19、nA = 0.8415 0.9093 0.1411 -0.7568 -0.9589 -0.2794 0.6570 0.9894 0.4121A_sinM = -0.6928 -0.2306 0.2316 -0.1724 -0.1434 -0.11430.3479 -0.0561 -0.4602 5.3 奇异值分解5.3.1 奇异值分解和矩阵结构5.3.1.1 奇异值分解定义5.3.1.2 矩阵结构的奇异值分解描述（1）秩 （2）零空间（Null space）和值空间（Range space）（3）范数（Norm）（4）矩阵条件数（5）子空间夹角（6）广义逆（Moore-Penrose pseu
20、doinverse）5.3.2 线性二乘问题的解5.3.2.1 矩阵除运算的广义化5.3.2.2 线性模型的最小二乘解【*例5.4.2.2-1】对于超定方程，进行三种解法比较。其中取MATLAB库中的特殊函数生成。（1）生成矩阵及，并用三种方法求解A=gallery(5);A(:,1)=;y=1.7 7.5 6.3 0.83 -0.082;x=inv(A*A)*A*y,xx=pinv(A)*y,xxx=Ay Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 7.087
21、751e-018.x = 3.4811e+000 5.1595e+000 9.5342e-001 -4.6562e-002xx = 3.4759e+000 5.1948e+000 7.1207e-001 -1.1007e-001Warning: Rank deficient, rank = 3 tol = 1.0829e-010.xxx = 3.4605e+000 5.2987e+000 0 -2.9742e-001 （2）计算三个解的范数nx=norm(x),nxx=norm(xx),nxxx=norm(xxx) nx = 6.2968e+000nxx = 6.2918e+000nxxx =
22、 6.3356e+000 （3）比较三种解法的方程误差e=norm(y-A*x),ee=norm(y-A*xx),eee=norm(y-A*xxx) e = 1.1020e+000ee = 4.7424e-002eee = 4.7424e-002 5.4 函数的数值导数和切平面5.4.1 法线【*例5.5.1-1】曲面法线演示。y=-1:0.1:1;x=2*cos(asin(y);%旋转曲面的“母线”X,Y,Z=cylinder(x,20);%形成旋转曲面surfnorm(X(:,11:21),Y(:,11:21),Z(:,11:21);%在曲面上画法线view(120,18)%控制观察角 图
23、 5.5.1-1 旋转半椭球面和法线5.4.2 偏导数和梯度5.4.2.1 理论定义5.4.2.2 数值计算指令【*例5.5.2.2-1】用一个简单矩阵表现diff和gradient指令计算方式。F=1,2,3;4,5,6;7,8,9Dx=diff(F)%相邻行差分Dx_2=diff(F,1,2)%相邻列差分。第三输入宗量2表示“列”差分。FX,FY=gradient(F)%数据点步长默认为1FX_2,FY_2=gradient(F,0.5)%数据点步长为0.5 F = 1 2 3 4 5 6 7 8 9Dx = 3 3 3 3 3 3Dx_2 = 1 1 1 1 1 1FX = 1 1 1
24、1 1 1 1 1 1FY = 3 3 3 3 3 3 3 3 3FX_2 = 2 2 2 2 2 2 2 2 2FY_2 = 6 6 6 6 6 6 6 6 6 【*例 5.5.2.2-2】研究偶极子(Dipole)的电势（Electric potential）和电场强度（Electric field density）。设在处有电荷，在处有电荷。那么在电荷所在平面上任何一点的电势和场强分别为，。其中。又设电荷，。clear;clf;q=2e-6;k=9e9;a=1.5;b=-1.5;x=-6:0.6:6;y=x;X,Y=meshgrid(x,y);%设置坐标网点rp=sqrt(X-a).2+
25、(Y-b).2);rm=sqrt(X+a).2+(Y+b).2);V=q*k*(1./rp-1./rm);%计算电势Ex,Ey=gradient(-V);%计算场强AE=sqrt(Ex.2+Ey.2);Ex=Ex./AE;Ey=Ey./AE;%场强归一化，使箭头等长cv=linspace(min(min(V),max(max(V),49);%产生49个电位值contourf(X,Y,V,cv,k-)%用黑实线画填色等位线图%axis(square)%在Notebook中，此指令不用title(fontname隶书fontsize22偶极子的场),hold onquiver(X,Y,Ex,Ey,
26、0.7)%第五输入宗量0.7使场强箭头长短适中。plot(a,b,wo,a,b,w+)%用白线画正电荷位置plot(-a,-b,wo,-a,-b,w-)%用白线画负电荷位置xlabel(x);ylabel(y),hold off 图 5.5.2.2-1 电偶极子的场和等位线5.5 函数的零点5.5.1 多项式的根5.5.2 一元函数的零点5.5.2.1 利用MATLAB作图指令获取初步近似解5.5.2.2 任意一元函数零点的精确解【*例 5.6.2.2-1】通过求的零点，综合叙述相关指令的用法。（1）构造一个内联函数对象被解函数以为自变量，和为参数。假如在fzero中直接采用字符串表示被解函数
27、，容易出错。因此先构造内联函数如下：y=inline(sin(t)2*exp(-a*t)-b*abs(t),t,a,b);% （2）作图法观察函数零点分布a=0.1;b=0.5;t=-10:0.01:10;%对自变量采样，采样步长不宜太大。y_char=vectorize(y);%为避免循环，把y改写成适合数组运算形式。 Y=feval(y_char,t,a,b);%在采样点上计算函数值。clf,plot(t,Y,r);hold on,plot(t,zeros(size(t),k);%画坐标横轴xlabel(t);ylabel(y(t),hold off 图 5.6.2.2-1 函数零点分布观
28、察图（3）利用zoom和ginput指令获得零点的初始近似值（在MATLAB指令窗中进行）zoom on%在MATLAB指令窗中运行，获局部放大图tt,yy=ginput(5);zoom off%在MATLAB指令窗中运行，用鼠标获5个零点猜测值。图 5.6.2.2-2 局部放大和利用鼠标取值图tt%显示所得零点初始猜测值（该指令可在Notebook中运行）。 tt = -2.0032 -0.5415 -0.00720.58761.6561 （4）求靠近tt(4) 的精确零点t4,y4,exitflag=fzero(y,tt(4),a,b)% Zero found in the interva
29、l: 0.57094, 0.60418.t4 =0。5993y4 = 0exitflag = 1 （5）求在tt(3)附近的精确零点从理论分析可知，是函数的一个零点。但即便是以十分靠近该零点的为搜索的初始值，也找不到，而却找到了另一个零点。原因是曲线没有穿越横轴。请看下面指令的运行结果。t3,y3,exitflag=fzero(y,tt(3),a,b) Zero found in the interval: 0.58266, -0.59706.t3 = -0.5198y3 = 0exitflag = 1 （6）观察fzero所采用的options缺省设置，并更改控制计算精度的相对误差设置。op
30、=optimset(fzero)%提取fzero所采用的options缺省设置 op = ActiveConstrTol: Display: final TolX: 2.2204e-016 TypicalX: op=optimset(tolx,0.01);%把终止计算的相对误差阈值设置得较大op.TolX%观察新设置值。注意TolX字母的大小写。 ans =0.0100 （7）利用新设置的选项参数重新求tt(4)附近的零点，以便比较。t4n,y4n,exitflag=fzero(y,tt(4),op,a,b)%采用新的op设置参数。 Zero found in the interval: 0.
31、57094, 0.60418.t4n =0。6042y4n =0。0017exitflag = 1 说明l 本例是采用内联函数形式求取函数零点的。l 若采用如下M函数文件（该文件必须放在搜索路径上）y_M.mfunction y=y_M(t,a,b)y=sin(t).2.*exp(-a*t)-b*abs(t);则相应的求零点指令是t4m,y4m,exitflag=fzero(y_M,tt(4),a,b)l 若直接用字符串表达函数，则应把被解函数自变量t改成x，参数a、b改成P1、P2。相应的具体指令如下P1=0.1;P2=0.5;y_C=sin(x).2.*exp(-P1*x)-P2*abs(
32、x);t4c,y4c,exitflag=fzero(y_C,tt(4),P1,P2)5.5.3 多元函数的零点【例 5.6.3-1】求解二元函数方程组的零点。（0）从三维坐标初步观察两函数图形相交情况x=-2:0.05:2;y=x;X,Y=meshgrid(x,y);%产生x-y平面上网点坐标F1=sin(X-Y);F2=cos(X+Y);F0=zeros(size(X);surf(X,Y,F1),xlabel(x),ylabel(y),view(-31,62),hold on,surf(X,Y,F2),surf(X,Y,F0),shading interp,hold off 图 5.6.3-
33、0 两函数的三维相交图（1）在某区域观察两函数0等位线的交点情况clear;x=-2:0.5:2;y=x;X,Y=meshgrid(x,y);%产生x-y平面上网点坐标F1=sin(X-Y);F2=cos(X+Y);v=-0.2, 0, 0.2;%指定三个等位值，是为了更可靠地判断0等位线的存在。contour(X,Y,F1,v)%画F1的三条等位线。hold on,contour(X,Y,F2,v),hold off%画F2的三条等位线。 图 5.6.3-1 两个二元函数0等位线的交点图（2）从图形获取零点的初始近似值在图5.6.3-1中，用ginput获取两个函数0等位线（即三线组中间那条
34、线）交点的坐标。x0,y0=ginput(2); %在图上取两个点的坐标disp(x0,y0) -0.7926 -0.7843 0.7926 0.7843 （3）利用fsolve求精确解。以求（0.7926,7843）附近的解为例。本例直接用字符串表达被解函数。注意：在此，自变量必须写成x(1), x(2)。假如写成xy(1), xy(2)，指令运行将出错。fun=sin(x(1)-x(2),cos(x(1)+x(2);%xy=fsolve(fun,x0(2),y0(2)% xy = 0.7854 0.7854 （4）检验fxy1=sin(xy(1)-xy(2);fxy2=cos(xy(1)+
35、xy(2);disp(fxy1,fxy2) 1.0e-006 * -0.0994 0.2019 说明l 指令可用以下任何一组指令取代。（A）内联函数形式指令fun=inline(sin(x(1)-x(2), cos(x(1)+x(2), x);%项x必须有。xy=fsolve(fun,x0(2), y0(2);（B）M函数文件形式及指令先用如下fun.m表示被解函数（并在搜索路径上）FUN.Mfunction ff=fun(x)ff(1)=sin(x(1)-x(2);ff(2)=cos(x(1)+x(2);然后运行指令xy=fsolve(fun,x0(2),y0(2) 。l 第四步检验中的结果
36、表明：所找零点处的函数值小于，是一个十分接近零的小数。该精度由options.TolFun控制。 options.TolFun的缺省值是1.0000e-006。它可以用下列指令看到options=optimset(fsolve);options.TolFun ans = 1.0000e-006 5.6 函数极值点5.6.1 一元函数的极小值点5.6.2 多元函数的极小值点【*例5.7.2-1】求的极小值点。它即是著名的Rosenbrocks Banana 测试函数。该测试函数有一片浅谷，许多算法难以越过此谷。（0）从三维等位线图初步观察测试函数x=-3:0.1:3;y=-2:0.1:4;X,Y
37、=meshgrid(x,y);F=100*(Y-X.2).2+(1-X).2;contour3(X,Y,F,300),xlabel(x),ylabel(y),axis(-3,3,-2,4,0,inf),view(161,22)hold on,plot3(1,1,0,.r,MarkerSize,20),hold off 图 5.7.2-1-0 三维等位线图显示出一条香蕉形浅谷（1）本例采用内联函数表示测试函数如下ff=inline(100*(x(2)-x(1)2)2+(1-x(1)2,x); （2）用单纯形法求极小值点x0=-1.2,1;sx,sfval,sexit,soutput=fminse
38、arch(ff,x0) Optimization terminated successfully: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-004 sx = 1.0000 1.0000sfval = 8.1777e-010sexit = 1soutput = iterations: 85 funcCount: 159 algorithm: Nelder-Mead simplex direct search （3）用拟牛顿法求极小值点ux,sfval,uexit,uoutput,grid,hess=fminunc(ff,x0) Warning: Gradient must be provided for trust-region method; using line-search method instead.