第四章数值计算

《第四章数值计算》由会员分享，可在线阅读，更多相关《第四章数值计算（31页珍藏版）》请在金锄头文库上搜索。

1、第四章第四章 数值计算数值计算 主要内容：主要内容：4.1 数值微积分数值微积分4.2 矩阵和代数方程矩阵和代数方程4.4 多项式运算多项式运算4.1数值微积分数值微积分4.1.1 近似数值极限及导数近似数值极限及导数 在在MATLAB数值计算中，既没有专门的求极限指令，数值计算中，既没有专门的求极限指令，也没有专门的求导指令。但也没有专门的求导指令。但MATLAB提供了与提供了与“求求导导”概念有关的概念有关的“求差分求差分”指令。指令。ndx=diff(X) 计算向量计算向量X的前向差分的前向差分ndx=diff(X,N) 计算向量计算向量X的的N阶前向差分阶前向差分nFX=gradien

2、t(F) 求一元（函数）梯度求一元（函数）梯度例例 4.1-3clfd=pi/100;t=0:d:2*pi;x=sin(t);dxdt_diff=diff(x)/d;dxdt_grad=gradient(x)/d;subplot(1,2,1)plot(t,x,b)hold onplot(t,dxdt_grad,m,LineWidth,8)plot(t(1:end-1),dxdt_diff,.k,MarkerSize,8)axis(0,2*pi,-1.1,1.1)title(0, 2pi)已知已知，采用，采用 diff 和和 gradient 计算算该函数在区函数在区间 中的近似导函数。中的近似

3、导函数。legend(x(t),dxdt_grad,dxdt_diff,Location,North)xlabel(t),box offhold offsubplot(1,2,2)kk=(length(t)-10):length(t);hold onplot(t(kk),dxdt_grad(kk),om,MarkerSize,8)plot(t(kk-1),dxdt_diff(kk-1),.k,MarkerSize,8)title(end-10, end)legend(dxdt_grad,dxdt_diff,Location,SouthEast)xlabel(t),box offhold off

4、 注意：不要企图借助数值计算求取极限注意：不要企图借助数值计算求取极限4.1.2 数值求和与近似数值积分数值求和与近似数值积分nsx=sum(X) 沿列方向求和沿列方向求和nscs=cumsum(X) 沿列方向累计求和沿列方向累计求和nst=trapz(x,y) 采样梯形法沿列方向求函数采样梯形法沿列方向求函数y 关于自变量关于自变量x的积分的积分nsct=cumtrapz(x,y) 采样梯形法沿列方向求函数采样梯形法沿列方向求函数 y 关于自变量关于自变量x的累计积分的累计积分clear d=pi/8; t=0:d:pi/2; y=0.2+sin(t); s=sum(y); s_sa=d*s

5、; s_ta=trapz(t,y); %图形显示指令图形显示指令图图 4.1-4 sum 和和trapz求积模式示意求积模式示意 说明：说明： （1）本例中的指令）本例中的指令s_ta=d*trapz(y); 可以用可以用s_ta=trapz(t,y)替换，即计算由替换，即计算由t，y所绘折线下的面积。所绘折线下的面积。（2）s=sum(y); s_sa=d*s; 用作近似积分是错误的。阶梯用作近似积分是错误的。阶梯虚线所占的自变量区间比积分区间多一个采样子区间。虚线所占的自变量区间比积分区间多一个采样子区间。即即不能把不能把d *sum(y)看作看作“矩形近似积分矩形近似积分”。（3）实际使

6、用中，应该把子区间划分得相当小以便获取较）实际使用中，应该把子区间划分得相当小以便获取较高精度的近似积分。显然，采样点愈多，积分精度愈高，高精度的近似积分。显然，采样点愈多，积分精度愈高，但精度无法定量确定。但精度无法定量确定。 4.1.3 计算精度可控的数值积分计算精度可控的数值积分 数值积分有闭型（数值积分有闭型（closedtype）算法、开型）算法、开型（opentype）算法，其主要区别在于：是否需）算法，其主要区别在于：是否需要计算积分区间端点处的函数值。要计算积分区间端点处的函数值。nS1=quad(fun,a,b,tol) Simpson法法nS1=quadl(fun,a,b,

7、tol) Lobatto法法nS2=dblquad(fun,xmin,xmax,ymin,ymax,tol)nS3=triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol)说明：说明：（1）fun是被积函数，可以是是被积函数，可以是字符串字符串、内联对象内联对象、匿匿名函数名函数和和M函数文件的函数句柄函数文件的函数句柄。（2）要保证对于向量形式的自变量（）要保证对于向量形式的自变量（一般为字母一般为字母x）输入，输出为长度相同的函数值向量。输入，输出为长度相同的函数值向量。（3）a，b为一重积分的下限和上限；为一重积分的下限和上限；xmin,xmax

8、,ymin,ymax,zmin,zmax是多重积分由是多重积分由内向外的积分限。内向外的积分限。（4）tol是个标量，用来控制绝对误差。默认时，积是个标量，用来控制绝对误差。默认时，积分的绝对精度为分的绝对精度为10-6.（1e-6）注：在编写被积函数时要遵循注：在编写被积函数时要遵循“数组运算数组运算”规则。规则。syms xIsym=vpa(int(exp(-x2),x,0,1) Isym =.74682413281242702539946743613185 format longd=0.001;x=0:d:1;Itrapz=trapz(x,exp(-x.*x) Itrapz = 0.74

9、682407149919 fx=exp(-x.2);Ic=quad(fx,0,1,1e-8) Ic = 0.74682413285445 format longfx1=x.2 %采用字符串表示被积函数采用字符串表示被积函数fx2=inline(x.y) %采用内联对象表示被积函数采用内联对象表示被积函数fx3=(x,y)x.y %采用匿名函数表示被积函数采用匿名函数表示被积函数s1=dblquad(fx1,0,1,1,2)s2=dblquad(fx2,0,1,1,2)s3=dblquad(fx3,0,1,1,2)s4=dblquad(fx4,0,1,1,2)sn = 0.40546626724

10、351 function z=fx4(x,y)z=x.y4.1.4 函数极值的数值求解函数极值的数值求解 nx,fval,exitflag,output=fminbnd(fun,x1,x2,options) 求一元函数在区间（求一元函数在区间（x1,x2）中的极小值。）中的极小值。nx,fval,exitflag,output=fminsearch(fun,x0,options) 单纯形法求多元函数的极值点单纯形法求多元函数的极值点fun字符串、内联对象、匿名函数、函数文件句柄；字符串、内联对象、匿名函数、函数文件句柄；输出量输出量x,fval分别是极值点和相应的目标函数极值；分别是极值点和相

11、应的目标函数极值；输出量输出量exitflag若给出大于若给出大于0的数，说明成功搜索；的数，说明成功搜索；输出量输出量output给出具体的优化算法和迭代次数。给出具体的优化算法和迭代次数。【例例4.1-7】已知，已知， 在区间在区间 ，求函数的极小值。，求函数的极小值。三种方法：三种方法：（1）教科书方法（自学）（）教科书方法（自学）（2）优化算法（）优化算法（fminbnd指令）（指令）（3）图形法求极小值（）图形法求极小值（ginput指令）指令）方法一：方法一：yx=(x+pi)*exp(abs(sin(x+pi);xn0,fval,exitflag,output=fminbnd(y

12、x, -pi/2, pi/2) 方法二：方法二：xx=-pi/2:pi/200:pi/2;yxx=(xx+pi).*exp(abs(sin(xx+pi);plot(xx,yxx) % fplot(yx,-pi/2,pi/2)gridxx,yy=ginput(1) 【例例】求二元函数在原点附近的最大求二元函数在原点附近的最大值：注：注：fminsearch被优化函数被优化函数fun中的多元自变量应采用单一变量名中的多元自变量应采用单一变量名的向量形式表达。的向量形式表达。ff=(x)(x(1)4+x(2)4-4*x(1)*x(2)-5); x0=0 0; x,g=fminsearch(ff,x0

13、) nt,Y=ode45(odefun,tspan,y0)采用四阶采用四阶RungeKutta数值积分法解微分方程数值积分法解微分方程注：注：odefun是待解微分方程的函数文件句柄。该函数文件的是待解微分方程的函数文件句柄。该函数文件的输出必须是待解函数的一阶导数。不管原问题是不是一阶输出必须是待解函数的一阶导数。不管原问题是不是一阶微分方程组，当使用微分方程组，当使用ode45求解时，必须转化成一阶微分求解时，必须转化成一阶微分方程组形式。方程组形式。式中，式中，y是（是（n*1）向量。）向量。tspan 定义求数值解的时间区间；定义求数值解的时间区间；y0 是一阶微分方程组的（是一阶微分

14、方程组的（n*1）初值列向量；）初值列向量；4.1.5 常微分方程的数值解常微分方程的数值解注意：当求解一阶微分注意：当求解一阶微分方程时，方程时，odefun可以是可以是待解微分方程的内联对待解微分方程的内联对象，匿名函数或函数文象，匿名函数或函数文件句柄。件句柄。【例例】解微分方程解微分方程odefun=inline(y-2*t/y,t,y);t,y=ode45(odefun,0,4,1);plot(t,y)% ode45(odefun,0,4,1)【例例4.1-9】求微分方程求微分方程在初始条件在初始条件情况下的解，并情况下的解，并图示。示。 （1）把高阶微分方程改写成一阶微分方程）把高

15、阶微分方程改写成一阶微分方程令令 ，则，则（2）根据上述方程组，编写）根据上述方程组，编写M函数文件函数文件dydt.mfunction ydot=dydt(t,y)ydot=y(2);2*(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)xlabel(t),title(x(t) % ode45(dydt,tspan,y0)4.2 矩阵和代数方程矩阵和代数方程4.2.1 矩阵运算和特征参数矩阵运算和特征参数一、矩阵运算一、矩阵运算2.转置运算转置运算【例

16、例 4.2-2】nA1=A 共轭转置共轭转置nA2=A. 非共轭转置非共轭转置二二 、矩阵的标量特征参数、矩阵的标量特征参数 MATLAB用来计算（限于大学教材中涉及的）矩用来计算（限于大学教材中涉及的）矩阵特征参数的指令，见表阵特征参数的指令，见表4.2-2.(P141)【例例4.2-3】矩阵标量特征参数计算示例。矩阵标量特征参数计算示例。【例例4.2-4】矩阵标量特征参数的性质。矩阵标量特征参数的性质。 4.2.2 矩阵变换和特征值分解矩阵变换和特征值分解nR,ci=rref(A) 借助初等变换把借助初等变换把A变换成行阶梯矩阵变换成行阶梯矩阵nX=null(A) A矩阵零空间的全部正交基

17、，满足矩阵零空间的全部正交基，满足AX0。 （求通解（齐次解）（求通解（齐次解）nZ=orth(A) A矩阵值空间的全部正交基，满足矩阵值空间的全部正交基，满足 span（Z）span（A）nV,D=eig(A) A矩阵的特征值、特征向量分解，使矩阵的特征值、特征向量分解，使 AV VD4.2.3 线性方程的解线性方程的解一、线性方程解的一般结论一、线性方程解的一般结论 对于含对于含n个未知数的个未知数的m个方程构成的方程组个方程构成的方程组AmnXb，它的解有一下几种可能：，它的解有一下几种可能：n1.当向量当向量b在矩阵在矩阵A列向量所张空间中，有准确解。列向量所张空间中，有准确解。ran

18、k(A,b)=rank(A)（1）若）若n=r，则解唯一。（，则解唯一。（ r= rank(A) ）（2）若）若nr，则解不唯一。（，则解不唯一。（ r= rank(A) ）n2.当向量当向量b不在矩阵不在矩阵A列向量所张空间中，则无列向量所张空间中，则无准确解，但存在最小二乘解。准确解，但存在最小二乘解。（了解）（了解）（1）当）当A列满秩时，存在唯一的最小二乘解列满秩时，存在唯一的最小二乘解（2）当）当A列不满秩时，存在最小范最小二乘解列不满秩时，存在最小范最小二乘解和最少非零元素最小二乘解。和最少非零元素最小二乘解。二、除法运算解方程二、除法运算解方程nAX=B,则则x=AB %左除（左

19、除（A-1*B）nXC=D,则则x=D/C %右除（右除（D*C-1）三、矩阵逆三、矩阵逆 A_1=inv(A) 求非奇异方阵求非奇异方阵A的逆。的逆。u计算结果表明：除法求解不但速度快，计算结果表明：除法求解不但速度快，而且精度高得多。而且精度高得多。【例例4.2-8】求方程求方程 的解。的解。 4.2.4 一般代数方程的解一般代数方程的解（1）利用）利用MATLAB作图指令（作图指令（ginput）获取初步近似解。）获取初步近似解。（2）利用）利用MATLAB指令求精确解。指令求精确解。x,favl=fzero(fun,x0) 求一元函数零点指令的最简格式求一元函数零点指令的最简格式x,f

20、avl=fsolve(fun,x0) 解非线性方程组的最简单格式解非线性方程组的最简单格式fun字符串、内联对象、匿名函数、函数文件句柄。被解函字符串、内联对象、匿名函数、函数文件句柄。被解函数的自变量一般采用字母数的自变量一般采用字母x。x0-表示零点初始猜测。表示零点初始猜测。【例例】求函数求函数 在（在（-2，-0.1）内的零点。）内的零点。fun=x*sin(x2-x-1);fzero(fun,-2,-0.1)?fplot(fun,-2,-0.1)grid onfzero(fun,-1,-1.2) %fzero(fun,-1.6)fzero(fun,-1.2,-0.1) %fzero(fun,-0.6)【例例】使用使用fsolve函数求解非线性方程组的解。函数求解非线性方程组的解。function f=fxx(x)a=x(1); b=x(2); c=x(3);f(1)=a2+b+sin(c);f(2)=a*b+c;f(3)=cos(a)+b2+2*c;x0=1,1,1;fsolve(fxx,x0)