MATLAB第4章.ppt

第第4章章 数值计算数值计算1本章学习目标n掌握LU分解和恰定方程组n掌握矩阵特征值和特征向量n掌握函数的极值点和数值积分的求解n理解随机数据的生成函数n理解傅立叶变换的各种性质，熟悉基本信号的频域转换，掌握应用FFT对典型信号进行频谱分析的方法n掌握龙格库塔法求解常微分方程数值解的方法。2主要内容n4.1 LU分解和恰定方程组的解n4.2 矩阵特征值和特征向量n4.3 函数的零点n4.4 函数的极值点n4.5 数值积分n4.6 随机数据的统计描述n4.7 傅里叶分析n4.8 常微分方程34.1 LU分解和恰定方程组的解主要调用格式：nL，U=lu(X)，其中是任意方阵，L是下三角矩阵，U是上三角矩阵，满足的条件式为；nL，U，P=lu(X)，其中是任意方阵，L是下三角矩阵，U是上三角矩阵，是置换矩阵，满足的条件式为PX=LU。n=lu(X)，其中是任意方阵，把下三角矩阵和上三角矩阵合并在矩阵中给出，满足的条件式为Y=L+U-I，该命令将损失置换矩阵的信息。注意：进行lu分解时，矩阵X必须是方阵。【例4-1】已知矩阵A=5-1 2;1 2-1;3 1 4，求其LU分解。A=5-1 2;1 2-1;3 1 4;L，U=lu(A)L*U%测试分解的正确性【例4-2】已知方程组 ，应用LU分解法求方程组的解。A=5-1 2;1 2-1;3 1 4;L，U=lu(A)b=2 3 5 x=U(Lb)4.1.2恰定方程组的解恰定方程组的解对于方程Ax=b，A为nm矩阵，有三种情况：n当n=m时，此方程成为“恰定”方程；n当nm时，此方程成为“超定”方程；n当n A=5-1 2;1 2-1;3 1 4;b=2 3 5;x=inv(A)*b方法二：左除法 A=5-1 2;1 2-1;3 1 4;b=2 3 5;x=Ab4.2 矩阵特征值和特征向量主要调用格式有：n=eig(A)，求矩阵的全部特征值，构成向量并进行输出。nV，D=eig(A)，求矩阵的特征向量矩阵和特征对角阵D，满足等式=的对角元素为A的特征值，而矩阵D的每一列向量为其所对应的特征向量。nV，D=eig(，nobalance)，当矩阵A中有截断误差数量级相差不大时，该指令更加精确。nobalance起误差调节作用。【例4-4】计算阶pascal矩阵的特征值和特征向量。A=pascal(4);E=eig(A)V，D=eig(A)4.3函数的零点n函数的零点是使f(x)=0的实数x，用二维图形表示也就是函数y=f(x)的图形和x轴交点的横坐标，即y=f(x)中y=0时x的取值。对于任意求解函数，可能有一个解，也可能有多个解，因此没有一个通用的解法。一般来说，可以通过作图，观察零点所在的位置，然后通过计算不断地逼近真实值，并在一定精度时停止计算。4.3.1一元函数的零点一元函数的零点在MATLAB中求解一元函数的零点的函数是fzero()，其调用格式如下：nx=fzero(fun，x0)，参数fun为一元函数名，x0表示求解的初始数值。nx,fval,exitflag,output=fzero(fun,x0,options)，options是优化迭代所采用的参数选项，在输出参数中，fval表示对应的函数值，exitflag表示程序退出的类型，output为优化信息的输出变量。以上函数的功能是求函数的零点x和该点的函数值fval，其中fun是待求解的一元函数，可以是字符串、内联函数、M函数文件名的函数句柄。给出初始值x0或区间位置，反复计算逼近真实值。fzero不仅能寻找零点，它还可以寻找函数等于任何常数值的点。例如，为了寻找f(x)=c的点，定义函数g(x)=f(x)-c，然后，在fzero中使用g(x)，就会找出g(x)为零的x值，它发生在f(x)=c时。【例4-5】求函数f(x)=x2-2x-5的零点。第一步，绘制函数的图形。%计算函数数值x=-3:0.1:5;y=x.2-2.*x-5;%绘制函数图形 plot(x，y);%锁定当前图形并添加网格线 hold on;grid on%绘制水平线 line(-3 5，0 0);%添加坐标名称xlabel(x);ylabel(y(x);第二步，利用ginput函数获得零点的初始近似值 xx，yy=ginput(2)%在MATLAB指令窗中运行，用鼠标获2个零点猜测值。第三步，计算xx(1)和xx(2)的精确零点 y=inline(x.2-2.*x-5，x);构造内联函数求解 x1，y1=fzero(y，xx(1)x2，y2=fzero(y，xx(2)说明：n上面是通过构造内联函数求解的，fzero()中被解函数还可以采用字符串或M函数文件名的函数句柄来进行求解。nfzero()函数只能返回一个局部零点，不能寻找到所有的零点。nfzero()函数的收敛速度与开始点的选取或区间的选取有关，一般是先找到零点的大概范围，然后在该范围内寻找零点。4.3.2 多元函数的零点多元函数的零点多元函数的零点比求解一元函数的零点要复杂，通常含有很多解，只能通过计算得到在某个值附件最接近的解。对于形式比较简单的二元函数，可以采用类似一元函数的方法，首先作图观察零点的范围，再通过函数计算真实值。对于复杂的多元函数，一般采用求解非线性方程组的方法求解。求解多元函数的命令是fsolve()，具体调用格式如下：nx=fsolve(fun，x0)nx，fval，exitflag，output=fsolve(fun，x0，options)其参数与fzero()函数的参数相似。例4-6求函数f(x)=x2-2y的零点。第一步，绘制函数的图形，观察零点的大概位置。%生成三维图形的数据网格x=-5:0.1:5;X，Y=meshgrid(x);%指定三个等位值，是为了更可靠地判断0等位线的存在。v=-0.3，0，0.3;%绘制函数的三条等高线图。contour(X，Y，Z，v)第二步，从图形获取一个零点的初始近似值用ginput()函数获取0等高线上一个接近零点的坐标。如图4-2 所示，在中间等高线上的点都是零点解。x，y=ginput(1)图4-2多元函数等高线图 第三步，利用fsolve求精确解。%编写函数文件，文件名为zero1.m。function f=zero1(x)f=x(1)2-2*x(2)%计算零点值 xout=fsolve(zero1，x，y)xout=1.5843%计算f(x)的值 f=xout(1)2-2*xout(2)计算结果表明，所找零点处的函数值小于，是一个十分接近零的小数。该精度由options.TolFun控制。options.TolFun的缺省值是1.0000e-006。它可以用下列指令看到 options=optimset(fsolve);options.TolFunans=1.0000e-0064.4函数的极值点n求函数的极值点是工程实践中常见的问题。因为求f(x)的极小值等价于求-f(x)的极大值，所以本节主要讨论极小值的问题。4.4.1 一元函数的极小值点一元函数的极小值点fminbnd()函数可以用来进行有约束的一元函数最小值的求解，调用格式如下：nx=fminbnd(fun，x1，x2)，求解函数在区间x1 x2的最小值，fun是目标函数，其表达方式同fzero和fsolve函数。nx=fminbnd(fun，x1，x2，options)，options为优化参数选项，由函数optimset()设定。nx，fval=fminbnd()，fval为目标函数的最小值。nx，fval，exitflag=fminbnd()，exitflag为终止迭代的原因。nx，fval，exitflag，output=fminbnd()，output为优化信息inbnd 进行有约束的一元函数最小值。例4-7求f(x)=x3-5x2+2x-5的极小值。%绘制函数曲线x=0:0.1:6;f=x.3-5*x.2+2*x-5;plot(x，f);%估计极小值点的坐标。由图4-3可以看出，在3 4间有极小值 xx，ff=ginput(1)xx=3.1382ff=-17.1898%编写函数文件，保存文件名为myfun1.mfunction f=myfun1(x)f=x.3-5*x.2+2*x-5;%求极小值 x，fval=fminbnd(myfun1，3，4)图4-3 一元函数的极小值4.4.2 多元函数的极小值点多元函数的极小值点多元函数的极小值点的解法分为无约束条件和有约束条件2种。在MATLAB中，无约束条件求极小值点分为单纯形法和拟牛顿法两种方法，调用格式如下：nx=fminsearch(fun，x0)；%对函数fun，从x0开始搜索最优值，返回最优变量x；nx=fminsearch(fun，x0，options)；%可以用optimset设置一些优化选项；nx，fval=fminsearch(.)；%返回f为最优函数值；nx，fval，exitflag=fminsearch(.)；%exitflag与fminbnd类似；nx，fval，exitflag，output=fminsearch(.)；%output与fminbnd类似；x=fminunc(fun，x0，options)x，fval=fminunc(.)x，fval，exitflag=fminunc(.)x，fval，exitflag，output=fminunc(.)x，fval，exitflag，output，grad=fminunc(.)；%grad为解x处的梯度值；x，fval，exitflag，output，grad，hessian=fminunc(.)；%目标函数在解x处的海赛（Hessian）值。例4-8求%编写函数文件，保存文件名为myfun2.m function f=myfun2(x)f=x(1).2+2*x(2).2;%x和y分别用x(1)和x(2)表示 x0=1 2;%给出估计点x0%单纯形法求解 x1，fval1=fminsearch(myfun2，x0)x1=1.0e-004*0.2685 -0.1362fval1=1.0917e-009%拟牛顿法求解 x2，fval2=fminunc(myfun2，x0)Warning:Gradient must be provided for trust-region method;using line-search method instead.In fminunc at 241Optimization terminated:relative infinity-norm of gradient less than options.TolFun.x2=1.0e-007*-0.4396 0.0840fval2=2.0737e-015的极小值点。2 有约束条件有约束条件的多元函数的极小值可以通过fmincon()求解，调用格式如下：x=fmincon(fun，x0，A，b)x=fmincon(fun，x0，A，b，Aeq，beq)x=fmincon(fun，x0，A，b，Aeq，beq，lb，ub)x=fmincon(fun，x0，A，b，Aeq，beq，lb，ub，nonlcon)x=fmincon(fun，x0，A，b，Aeq，beq，lb，ub，nonlcon，options)例4-9求 的极值点，要求 。%编写函数文件mycon.mfunction f=mycon(x)f=x(1)2+2*x(2)2;%变换约束条件为，得到A和b，满足A=1，3;-1，-3;b=5;0;%赋初值x0=0.5;0.5%求极值和该点函数值x，fval=fmincon(mycon，x0，A，b)。4.5数值积分n在MATLAB中，可以使用多种方法实现微积分运算，如数值积分、符号积分、样条积分和Simulink积分器等。本节主要介绍数值积分。4.5.1 一元函数的数值积分一元函数的数值积分n在MATLAB中，一元函数的数值积分的函数有quad()和quadl()。quad()采用低阶的自适应Simpson法，quadl()采用高阶的自适应法Lobatto法。自适应Simpson法quad()采用低阶的自适应Simpson法，精度较高，经常使用，其调用格式如下：nq=quad(fun，a，b)nq=quad(fun，a，b，tol)nq=quad(fun，a，b，tol，trace)q=quad(fun，a，b，tol，trace，p1，p2)式中，fun是被积函数。a，b为积分的下限和上限。tol为积分允许误差。trace设置为1则显示递推过程的数据表，为0则不显示。p1，p2用来对函数fun直接传递参数。自适应法Lobatto法quadl()采用高阶的自适应法Lobatto法，精度更高，最常使用，其调用格式如下：nq=qual(fun，a，b)nq=quadl(fun，a，b，tol)nq=quad(fun，a，b，tol，trace)nq=quad(fun，a，b，tol，trace，p1，p2)例4-10求 。%使用内联函数inline fun=inline(1./(x.3-2*x-5),x);%自适应Simpson法计算积分 Q1=quad(fun，0，2)Q 1=-0.4605%自适应法Lobatto法计算积分 Q2=quad(F，0，2)Q2=-0.4605。4.5.2 多重数值积分多重数值积分二重数值积分函数的调用格式如下：nq=dblquad(fun，xmin，xmax，ymin，ymax)nq=dblquad(fun，xmin，xmax，ymin，ymax，tol)nq=dblquad(fun，xmin，xmax，ymin，ymax，tol，method)功能：调用函数在区域xmin，xmax，ymin，ymax上计算二元函数f(x，y)的二重积分。fun为被积函数。tol为允许误差，缺省时允许误差为10-6。Method为积分方法选项，缺省时为quad。例4-11计算方法一：编写匿名函数计算二重积分 Q=dblquad(x，y)(y*sin(x)+x*cos(y)，pi，2*pi，0，pi)Q=-9.8696方法二：编写函数文件计算二重积分%建立一个函数文件integrnd.mfunction z=integrnd(x，y)z=y*sin(x)+x*cos(y);%调用dblquad()函数求解 Q=dblquad(integrnd，pi，2*pi，0，pi)Q=-9.8696方法三：建立内联函数求解 fun=inline(y*sin(x)+x*cos(y)fun=Inline function:fun(x，y)=y*sin(x)+x*cos(y)Q=dblquad(fun，pi，2*pi，0，pi)Q=-9.86964.6随机数据的统计描述MATLAB有两个产生随机数据的命令。rand(m，n)产生在(0，1)之间均匀分布的m*n的随机数矩阵，其均值为0.5，标准差为0.2887;randn(m，n)产生服从正态分布的m*n的随机数矩阵，其均值为0，标准差为1。4.7 FOURIER分析n傅里叶变换在信号处理领域有及其重要的应用。它既可以对连续信号进行分析，也可以对离散信号进行分析。连续傅立叶变换实际上是计算傅立叶积分，将在第10章予以介绍，本节主要介绍离散傅里叶变换（Discrete Fourier Transform，缩写为DFT）。函数调用格式如下：nfft(X)，如果X是向量，返回其离散傅里叶变换，如果X是矩阵，函数对矩阵X的每一列进行离散傅里叶变换。nfft(X，N)，计算N点离散傅里叶变换。用N指定了傅里叶变换的长度，如果向量X的长度小于N，则在不足部分补0，使向量X的长度为N；若X的长度大于N，则删除超出N的那些元素。nfft(X，N，dim)，在dim维上进行傅里叶变换。当dim1时，该函数作用于X的每一行，当dim2时，该函数作用于X的每一列。表4-1傅里叶变换函数例4-12 求的离散傅里叶变换。N=1000，t在01s采样，绘制其时域及频域曲线。N=1000;%采样点数t=linspace(0，1，N);%时间序列 x=2*sin(2*pi*50*t)+0.5*sin(2*pi*180*t);%信号subplot(211);plot(1000*t(1:50)，x(1:50)%绘出时域曲线xlabel(时间);ylabel(样本);title(时域曲线);grid on;y=fft(x，1000);%对信号进行快速Fourier变换P=abs(y)*2/1000;%幅值f=1000*(0:500)/1000;%频率序列subplot(212);plot(f，P(1:501)%绘出随频率变化的振幅曲线xlabel(频率);ylabel(幅值);title(频域曲线);grid on;4.8常微分方程常微分方程函数调用格式为nt，y=solver(odefun，tspan，y0，options)nt，y，te，ye，ie=solver(odefun，tspan，y0，options)参数说明：solver为命令ode45、ode23，ode113，ode15s，ode23s，ode23t，ode23tb之一。odefun 为显式一阶常微分方程y=f(t，y)，若为高阶微分方程式，则必须将其化成一阶形式。tspan 为选择积分区间（即求解区间）的向量tspan=t0，tf。要获得问题在其他指定时间点t0，t1，t2，tf上的解，则tspan=t0，t1，t2，tf，其元素必须按单调次序排序。y0 为选择初始条件的向量。options 为设置算法参数，由命令odeset进行设置。t和y分布给出时间向量和相应的状态向量。参数te，ye，ie只有在设置”event”记录时，才有对应的输出数据。例4-13已知一阶微分方程 ，初始条件 ，求t=0:0.1:1时的y向量，并与解析解进行比较。第一步，建立函数文件funt.mfunction dy=funt(t，y)dy=y+2*t-1;第二步，求微分方程数值解%用ode45求微分方程数值解 t，y=ode45(funt，0:0.1:1，1);t%进行符号运算，求解析解 t=0:0.1:1;s=dsolve(Dy=y+2*t-1，y(0)=1，t)s=-2*t-1+2*exp(t)%求微分方程式的精确解 y1=-2*t-1+2*exp(t)%编写函数文件odet2.mfunction dy=odet2(t，y)dy=zeros(2，1);dy(1)=y(2);dy(2)=(1-y(1)2)*y(2)-y(1);%初始条件设定 tspan=0 20;y0=1 0;%求解常微分方程 t，y=ode45(odet2，tspan，y0)；%绘制计算结果并进行标注 plot(t，y)plot(t，y(:，1)，-，t，y(:，2)，-);title(常微分方程的解);legend(y1，y2);xlabel(t);ylabel(y);本章小结本章小结 数值计算在科学研究和工程中有十分广泛的应用，MATLAB的数值运算功能强大，并且计算方便。本章详细讲解了LU分解和恰定方程组的解、矩阵特征值和特征向量、函数的零点、函数的极值点、数值积分、随机数据的统计描述、FOURIER分析、常微分方程等内容。