第6章_MATLAB数值计算_part2

第6章MATLAB数值计算6.1数据处理与多项式计算6.2数值微积分6.3线性方程组求解6.4最优化问题求解6.5常微分方程的数值求解6.2数值微积分6.2.1数值微分（导数）不关心微分的形式和性质，只关心该微分在一串离散点的近似值以及所计算的近似值有多大的误差。MATLAB下求数值导数的两种方法：1）用多项式对任意函数进行拟合，再利用多项式求导函数进行求导，从而得出某点的导数值；2）用向前差商的方法求任意函数在某点的导数值。用向前差分方法求数值微分(1)在MATLAB中，没有直接提供求数值微分的函数，只有计算向前差分的函数diff，其调用格式为：DX=diff(X)：计算向量X的向前差分，DX(i)=X(i+1)-X(i)，i=1,2,…,n-1。DX=diff(X,n)：计算X的n阶向前差分。例如，diff(X,2)=diff(diff(X))。DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省状态)，按列计算差分；dim=2，按行计算差分。用向前差商方法求数值导数向前差分除以步长：diff(x)/h称为函数在所求的点处以h(h0)为步长的向前差商。只要步长足够小，向前差分即无限接近于函数在该点的微分，向前差商的结果无限接近于函数在该点的数值导数。还可以用向后差分、中心差分的方法求数值微分例6.18设x由[0,2π]间均匀分布的10个点组成，求sinx的1~3阶差分。命令如下：X=linspace(0,2*pi,10);Y=sin(X);DY=diff(Y);%计算Y的一阶差分D2Y=diff(Y,2);%计算Y的二阶差分，也可用命令diff(DY)计算D3Y=diff(Y,3);%计算Y的三阶差分，也可用diff(D2Y)或diff(DY,2)例6.19用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做出f'(x)的图像。程序如下：f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2');g=inline('(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2-x+12)/2+1/6./(x+5).^(5/6)+5');x=-3:0.01:3;p=polyfit(x,f(x),5);%用5次多项式p拟合f(x)dp=polyder(p);%对拟合多项式p求导数dpdpx=polyval(dp,x);%求dp在假设点的函数值dx=diff(f([x,3.01]))/0.01;%直接对f(x)求数值导数gx=g(x);%求函数f的导函数g在假设点的导数plot(x,dpx,x,dx,'.',x,gx,'-');%作图255122)(623xxxxxxf6.2.2数值积分数值积分方法求解定积分的数值方法多种多样，如简单的梯形法、辛普生(Simpson)法、牛顿－柯特斯(Newton-Cotes)法等都是经常采用的方法。基本思想将整个积分区间[a,b]分成n个子区间[xi,xi+1]，i=1,2,…,n，其中x1=a，xn+1=b；用区间的简单函数（插值多项式）替代被积函数，对每个区间求面积，并进行累加；求定积分问题分解为求和问题。公式：12()()()d()d()2()()d()d(()4()())62bbaabbaafafbfxxpxxbabaabfxxpxxfaffb11()()2()2nnkhTfafbfakh11/20111/201(()4()(1))6()4()2()()6nnkkkknnkkkkhSfxfxfxhfafxfxfb数值积分的实现1.被积函数是一个解析式quad(filename,a,b,tol,trace)quadl(filename,a,b,tol,trace)filename:被积函数名；a,b积分区间tol：积分精度；trace：控制是否展现积分过程例6.20用两种不同的方法求定积分。先建立一个函数文件ex.m：functionex=ex(x)ex=exp(-x.^2);然后在MATLAB命令窗口，输入命令：formatlongI=quad('ex',0,1)%注意函数名应加字符引号I=0.74682418072642I=quadl('ex',0,1)I=0.74682413398845也可不建立关于被积函数的函数文件，而使用语句函数(内联函数)求解，命令如下：g=inline('exp(-x.^2)');%定义一个语句函数g(x)=exp(-x^2)I=quadl(g,0,1)%注意函数名不加'号I=0.74682413398845formatshort102dxeIx2被积函数由一个表格定义在科学实验和工程应用中，函数关系往往是不知道的，只有实验测定的一组样本点和样本值，这时，就无法使用quad函数计算其定积分。在MATLAB中，对由表格形式定义的函数关系的求定积分问题用trapz(X,Y)函数。其中向量X、Y定义函数关系Y=f(X)。X、Y是两个等长的向量：X=(x1,x2,…,xn)，Y=(y1,y2,…,yn)，x1x2…xn积分区间是[x1，xn]。例6.21用trapz函数计算定积分。在MATLAB命令窗口，输入命令：X=0:0.01:1;Y=exp(-X.^2);trapz(X,Y)ans=0.7468102dxeIx3二重积分数值求解dblquad函数就可以直接求出二重定积分的数值解。该函数的调用格式为：I=dblquad(f,a,b,c,d,tol,trace)该函数求f(x,y)在[a,b]×[c,d]区域上的二重定积分。参数tol，trace的用法与函数quad完全相同。例6.22计算二重定积分。(1)建立一个函数文件fxy.m：functionf=fxy(x,y)globalki;ki=ki+1;%ki用于统计被积函数的调用次数f=exp(-x.^2/2).*sin(x.^2+y);(2)调用dblquad函数求解。globalki;ki=0;I=dblquad('fxy',-2,2,-1,1)kiI=1.57449318974494ki=1038112222/)sin(2dxdyyxeIx6.3离散傅立叶变换6.3.1离散傅立叶变换算法简要6.3.2离散傅立叶变换的实现一维离散傅立叶变换函数，其调用格式与功能为：(1)fft(X)：返回向量X的离散傅立叶变换。设X的长度(即元素个数)为N，若N为2的幂次，则为以2为基数的快速傅立叶变换，否则为运算速度很慢的非2幂次的算法。对于矩阵X，fft(X)应用于矩阵的每一列。(2)fft(X,N)：计算N点离散傅立叶变换。它限定向量的长度为N，若X的长度小于N，则不足部分补上零；若大于N，则删去超出N的那些元素。对于矩阵X，它同样应用于矩阵的每一列，只是限定了向量的长度为N。(3)fft(X,[],dim)或fft(X,N,dim)：这是对于矩阵而言的函数调用格式，前者的功能与FFT(X)基本相同，而后者则与FFT(X,N)基本相同。只是当参数dim=1时，该函数作用于X的每一列；当dim=2时，则作用于X的每一行。值得一提的是，当已知给出的样本数N0不是2的幂次时，可以取一个N使它大于N0且是2的幂次，然后利用函数格式fft(X,N)或fft(X,N,dim)便可进行快速傅立叶变换。这样，计算速度将大大加快。相应地，一维离散傅立叶逆变换函数是ifft。ifft(F)返回F的一维离散傅立叶逆变换；ifft(F,N)为N点逆变换；ifft(F,[],dim)或ifft(F,N,dim)则由N或dim确定逆变换的点数或操作方向。例6.23给定数学函数x(t)=12sin(2π×10t+π/4)+5cos(2π×40t)取N=128，试对t从0~1秒采样，用fft作快速傅立叶变换，绘制相应的振幅-频率图。在0~1秒时间范围内采样128点，从而可以确定采样周期和采样频率。由于离散傅立叶变换时的下标应是从0到N-1，故在实际应用时下标应该前移1。又考虑到对离散傅立叶变换来说，其振幅|F(k)|是关于N/2对称的，故只须使k从0到N/2即可。程序如下：N=128;%采样点数T=1;%采样时间终点t=linspace(0,T,N);%给出N个采样时间ti(I=1:N)x=12*sin(2*pi*10*t+pi/4)+5*cos(2*pi*40*t);%求各采样点样本值xdt=t(2)-t(1);%采样周期f=1/dt;%采样频率(Hz)X=fft(x);%计算x的快速傅立叶变换XF=X(1:N/2+1);%F(k)=X(k)(k=1:N/2+1)f=f*(0:N/2)/N;%使频率轴f从零开始plot(f,abs(F),'-*')%绘制振幅-频率图xlabel('Frequency');ylabel('|F(k)|')6.4线性方程组求解6.4.1直接解法1．利用左除运算符的直接解法对于线性方程组Ax=b，可以利用左除运算符“\”求解：x=A\b例6.24用直接解法求解下列线性方程组。命令如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';x=A\b04662975135243214324214321xxxxxxxxxxxxxx2．利用矩阵的分解求解线性方程组*矩阵分解是指根据一定的原理用某种算法将一个矩阵分解成若干个矩阵的乘积。常见的矩阵分解有LU分解、QR分解、Cholesky分解，以及Schur分解、Hessenberg分解、奇异分解等。分解后可以大大提高运算速度。例如：A=LU111211121222221200nnnnnnnnuuullluuLUlllu04662975135243214324214321xxxxxxxxxxxxxx例6.25用LU分解求解线性方程组。A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';[L,U]=lu(A);x=U\(L\b)[L,U,P]=lu(A);x=U\(L\P*b)(1)LU分解矩阵的LU分解就是将一个矩阵表示为一个交换下三角矩阵和一个上三角矩阵的乘积形式。线性代数中已经证明，只要方阵A是非奇异的，LU分解总是可以进行的。MATLAB提供的lu函数用于对矩阵进行LU分解，其调用格式为：[L,U]=lu(X)：使之满足X=LU。矩阵X必须是方阵[L,U,P]=lu(X)：U---上三角阵L---下三角阵P--置换矩阵PX=LU。X同样必须是方阵。线性方程组Ax=b的解x=U\(L\b)或x=U\(L\P*b)这样可以大大提高运算速度。例6.25用LU分解求解例6.24中的线性方程组。命令如下：A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];b=[13,-9,6,0]';[L,U]=lu(A);x=U\(L\b)或采用LU分解的第2种格式，命令如下：[L,U,P]=lu(A);x=U\(L\P*b)(2)QR分解对矩阵X进行QR分解，就是把X分解为一个正交矩阵Q和一个上三角矩阵R的乘积形式。QR分解只能对方阵进行。MATLAB的函数qr可用于对矩阵进行QR分解，其调用格式为：[Q,R]=qr(X)：X=QRQ正交矩阵R上三角矩阵[Q,R,E]=qr(X)：XE=QRQ正交矩阵、R上三角矩阵E置换矩阵实现QR分解后，线性方程组Ax=b的解x=R\(Q\b)或x=E(R\(Q\b))例6.26用QR分解求解例6.24中的线性方程组。命令如下：A=[2,1,-5,1;1