2第3-4讲-数值计算_New.doc

1、1第 4 章 数值计算与符号计算相比，数值计算在科研和工程中的应用更为广泛。本章主要讲授 MATLAB 的数值计算能力。4.1 数值微积分4.1.1 近似数值极限及导数【例 4.1-1】设 ， ，试用机器零阈值 eps 替代理论xfsin2co1)(xfsin)(0 计算极限 ， 。lim01fLxl02Lxx=eps;L1=(1-cos(2*x)/(x*sin(x),L2=sin(x)/x, L1 =2L2 =1 【例 4.1-2】已知 ，求该函数在区间 中的近似导函数。)sin(tx2 ,0d=pi/100;t=0:d:2*pi;x=sin(t);x_d=sin(t+d);dxdt=(x_

2、d-x)/d;plot(t,x,LineWidth,5)hold onplot(t,dxdt_d)hold offlegend(x(t),dx/dt)xlabel(t) 24.1.2 数值求和与近似数值积分【例 4.1-4 】求积分 ，其中 。 2/0)()(dtyxs )sin(2.0tycleard=pi/8;t=0:d:pi/2;y=0.2+sin(t);s=sum(y); %求和s_sa=d*s;s_ta=d*trapz(y); %梯形法求积分disp(sum 求得积分,blanks(3),trapz 求得积分)disp(s_sa, s_ta)sum 求得积分 trapz 求得积分1.

3、5762 1.3013图 4.1-4 sum 和 trapz 求积模式示意4.1.3 计算精度可控的数值积分【例 4.1-5 】求 。dxeI1023format long %15 位数字显示d=0.001;x=0:d:1;Itrapz=d*trapz(exp(-x.*x) %6 位数字有效Itrapz =0.74682407149919 fx=exp(-x.2); %字符函数Ic=quad(fx,0,1,1e-8) %指定精度 1e-8 计算结果达到 1e-10，Ic =0.74682413285445 【例 4.1-6 】求 。dxys2 10format longs_n=dblquad(

4、x,y)x.y,0,1,1,2) %二重积分，匿名函数，精度默认值 1e-6s_n =0.40546626724351 4.1.4 函数极值的数值求解【例 4.1-7】已知 ，在 区间，求函数的)sin()(xexy 2/x极小值。（1）x1=-pi/2;x2=pi/2;yx=(x)(x+pi)*exp(abs(sin(x+pi);xn0,fval,exitflag,output=fminbnd(yx,x1,x2) %求一元函数极小值（2）xx=-pi/2:pi/200:pi/2;yxx=(xx+pi).*exp(abs(sin(xx+pi);plot(xx,yxx)xlabel(x),gri

5、d on 4图 4.1-5 在-pi/2,pi/2区间中的函数曲线【例 4.1-8】求 的极小值点。它即是著名的22)1()(10),( xxyxf Rosenbrocks “Banana“ 测试函数，它的理论极小值是 。1,y（1）ff=(x)(100*(x(2)-x(1).2)2+(1-x(1)2); %求二元函数极小值（2）x0=-5,-2,2,5;-5,-2,2,5;sx,sfval,sexit,soutput=fminsearch(ff,x0) sx =1.0000 -0.6897 0.4151 8.08861.0000 -1.9168 4.9643 7.8004sfval =2.4

6、112e-010sexit =1soutput = iterations: 384funcCount: 615algorithm: Nelder-Mead simplex direct search %单纯形法message: 1x196 char （3）检查目标函数值format short edisp(ff(sx(:,1),ff(sx(:,2),ff(sx(:,3),ff(sx(:,4) 2.4112e-010 5.7525e+002 2.2967e+003 3.3211e+005 思考：如何求解极大值？4.1.5 常微分方程的数值解【例 4.1-9】求微分方程 ， ，在初始条件0)1(2

7、2xdtdtx25情况下的解，并图示。（见图 4.1-7 和 4.1-8）0)(,1)0(dtx（1）y1=x, y2=dx/dt, dy1/dt; dy2/dt=y2; mu*(1-y12)*y2-y1（2）DyDt.mfunction ydot=DyDt(t,y)mu=2;ydot=y(2);mu*(1-y(1)2)*y(2)-y(1);（3）解算微分方程tspan=0,30;y0=1;0;tt,yy=ode45(DyDt,tspan,y0);plot(tt,yy(:,1) %绘制 yy 第一列即 y1=x 的值xlabel(t),title(x(t) 图 4.1-7 微分方程解（4）pl

8、ot(yy(:,1),yy(:,2)xlabel(位移),ylabel( 速度) 图 4.1-8 平面相轨迹64.2 矩阵和代数方程4.2.1 矩阵运算和特征参数仅从数据排列看，矩阵就是二维数组。规定了特定的运算和意义的数组就成为矩阵。一 矩阵运算【例 4.2-2 】观察矩阵的转置操作和数组转置操作的差别。format ratA=magic(2)+j*pascal(2)A =1 + 1i 3 + 1i 4 + 1i 2 + 2i %A1=AA2=A. A1 =1 - 1i 4 - 1i 3 - 1i 2 - 2i A2 =1 + 1i 4 + 1i 3 + 1i 2 + 2i 二 矩阵的标量特

9、征参数【例 4.2-3】矩阵标量特征参数计算示例。A=reshape(1:9,3,3)r=rank(A)d3=det(A)d2=det(A(1:2,1:2)t=trace(A) A =1 4 7 2 5 8 3 6 9 r =2 d3 =0 d2 =-3 t =15 【例 4.2-4】矩阵标量特征参数的性质。rand(state,0)A=rand(3,3);B=rand(3,3);C=rand(3,4);D=rand(4,3);7tAB=trace(A*B) %矩阵位置交换秩不变tBA=trace(B*A) tCD=trace(C*D)tDC=trace(D*C) tAB =3.6697tBA

10、 =3.6697tCD =2.1544tDC =2.1544 d_A_B=det(A)*det(B)%同阶矩阵乘积行列式等于个矩阵行列式乘积dAB=det(A*B)dBA=det(B*A) d_A_B =0.0846dAB =0.0846dBA =0.0846 dCD=det(C*D)dDC=det(D*C) %非同阶矩阵乘积行列式随乘积顺序不同而不同dCD =-0.012362dDC =-2.1145e-018 4.2.2 矩阵的变换和特征值分解【例 4.2-7】简单实阵的特征值分解。 （1）A=1,-3;2,2/3V,D=eig(A) A =1.0000 -3.00002.0000 0.6

11、667V =0.7746 0.7746 0.0430 - 0.6310i 0.0430 + 0.6310iD =0.8333 + 2.4438i 0 0 0.8333 - 2.4438i （2）VR,DR=cdf2rdf(V,D) %把复数特征值对角阵 D 转换成实数块矩阵 DRVR =0.7746 00.0430 -0.63108DR =0.8333 2.4438-2.4438 0.8333 （3）A1=V*D/VA1_1=real(A1)A2=VR*DR/VRerr1=norm(A-A1,fro)err2=norm(A-A2,fro) A1 =1.0000 - 0.0000i -3.000

12、0 2.0000 + 0.0000i 0.6667 A1_1 =1.0000 -3.00002.0000 0.6667A2 =1.0000 -3.00002.0000 0.6667err1 =4.4409e-016err2 =4.4409e-016 4.2.3 线性方程的解一 线性方程解的一般结论二 除法运算解方程【例 4.2-8】求方程 的解。165432847306951x（1）A=reshape(1:12,4,3);b=(13:16); （2）ra=rank(A)rab=rank(A,b) %秩相等有准确解；A 的列数等于矩阵的秩有唯一解，否则解不唯一%秩不相等无准确解，有最小二乘解；%

13、列满秩存在唯一解，否则存在最小范数和最少非零元素最小二乘解ra =2rab =2 （3）9xs=Ab; %求特解xg=null(A); %求齐次方程解c=rand(1); %随机数ba=A*(xs+c*xg) %随机全解norm(ba-b) Warning: Rank deficient, rank = 2, tol = 1.8757e-014.ba =13.000014.000015.000016.0000ans =1.3874e-014 三 矩阵逆【例 4.2-9】“逆阵”法和“左除”法解恰定方程的性能对比（1）randn(state,0);A=gallery(randsvd,300,2e

14、13,2); %A 是条件数为 2e13 的 300 阶随机矩阵x=ones(300,1);b=A*x;cond(A) ans =2.0069e+013 （2）ticxi=inv(A)*b;ti=toceri=norm(x-xi)rei=norm(A*xi-b)/norm(b) ti =0.0341eri =0.0839rei =0.0048 （3）tic;xd=Ab;td=tocerd=norm(x-xd)red=norm(A*xd-b)/norm(b) td =100.0172erd =0.0127red =9.4744e-015 4.2.4 一般代数方程的解【例 4.2-10 】求 的零

15、点。 5.0)(sin)1.2tettft（1）数值计算y_C=inline(sin(t).2.*exp(-0.1*t)-0.5*abs(t),t); %内联对象函数t=-10:0.01:10;Y=y_C(t);clf,plot(t,Y,r);hold onplot(t,zeros(size(t),k);xlabel(t);ylabel(y(t)hold off 图 4.2-1 函数零点分布观察图zoom ontt,yy=ginput(5);zoom off 图 4.2-2 局部放大和利用鼠标取值图11tt tt =-2.0039-0.5184-0.00420.60521.6717 t4,y4

16、=fzero(y_C,0.1) %求一元函数零点t4 =0.5993y4 =1.1102e-016 t44,y44=fsolve(y_C,0.5) %解非线性方程组4.3 多项式运算和卷积4.3.1 多项式的运算函数一 多项式表达方式的约定二 多项式运算函数【例 4.3-1】求 的“商”及“余”多项式。1)(4)2(3s（1）format ratp1=conv(1,0,2,conv(1,4,1,1);p2=1 0 1 1;q,r=deconv(p1,p2);cq=商多项式为 ; cr=余多项式为 ;disp(cq,poly2str(q,s),disp(cr,poly2str(r,s)商多项式为

17、 s + 5余多项式为 5 s2 + 4 s + 3 【例 4.3-2 】矩阵和特征多项式，特征值和多项式根。（1）A=11 12 13;14 15 16;17 18 19;PA=poly(A) PPA=poly2str(PA,s) PA =1.0000 -45.0000 -18.0000 0.000012PPA =s3 - 45 s2 - 18 s + 2.163e-014 （2）s=eig(A)r=roots(PA) s =45.3965-0.39650.0000r =45.3965-0.39650.0000 【例 4.3-3 】构造指定特征根的多项式。R=-0.5,-0.3+0.4*i,

18、-0.3-0.4*i;P=poly(R)PR=real(P) PPR=poly2str(PR,x) P =1.0000 1.1000 0.5500 0.1250PR =1.0000 1.1000 0.5500 0.1250PPR =x3 + 1.1 x2 + 0.55 x + 0.125 【例 4.3-4 】多项式求值指令 polyval 与 polyvalm 的本质差别。（1）clearp=1,2,3;poly2str(p,x)X=1,2;3,4 ans =x2 + 2 x + 3X =1 23 4 （2）va=X.2+2*X+3Va=polyval(p,X) va =6 1118 27Va

19、 =6 1118 27 （3）vm=X2+2*X+3*eye(2)Vm=polyvalm(p,X) vm =12 141321 33Vm =12 1421 33 （4）cp=poly(X);poly2str(cp,x)cpXa=polyval(cp,X)cpX=polyvalm(cp,X) %哈密顿订立ans =x2 - 5 x - 2cpXa =-6 -8-8 -6cpX =1.0e-015 *0.2220 00 0.2220 4.3.2 多项式拟合和最小二乘法一 多项式拟合【例 4.4-5】 给定数据组 x0 , y0 ，求拟合三阶多项式，并图示拟合情况。（见图 4.4-1）x0=0:0.

20、1:1;y0=-.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22;n=3;P=polyfit(x0,y0,n)xx=0:0.01:1;yy=polyval(P,xx);plot(xx,yy,-b,x0,y0,.r,MarkerSize,20)legend(拟合曲线,原始数据,Location,SouthEast)xlabel(x) P =56.6915 -87.1174 40.0070 -0.904314图 4.4-1 采用三阶多项式所得的拟合曲线习题 41采用数值计算方法，画出 在 区间曲线，并计算 。dtxyx0sin)( 10

21、 , )5.4(y答案s45 =1.65410 2 4 6 8 1000.20.40.60.811.21.41.61.822用 quad 求取 的数值积分，并保证积分的绝对精度为 。dxexsin7.15 91015答案sq =1.08784993815498 3求函数 在区间 中的最小值点。5.0812cos5.)5(sin)206. ttettft ,答案最小值点是 -1.28498111480531相应目标值是 -0.186048010065454设 ，用数值法求 。0)(,1)0(,2)(3)(2 dtytydtty 5.0)(ty答案数值解y_05 =0.78958020790127 5求 的实数解。0sin10.2.tett答案t =2.7341yt =2.2204e-015 -1 0 1 2 3 4 5-12-10-8-6-4-20246现有一组实验数据 x, y（数据从 prob_data418.mat 获得），试求这组数据的 5 阶拟合多项式。答案P =-0.0039 0.0338 -0.0227 -0.4456 0.9590 -0.036416-1 0 1 2 3 4-1.4-1.2-1-0.8-0.6-0.4-0.200.20.40.6