matlab数值计算与符号计算.ppt

第4章MATLAB数值计算与符号计算4.1曲线拟合与插值运算4.2数值微积分4.3线性方程组求解4.4常微分方程的数值求解4.5MATLAB符号计算4.6级数4.7实验五数值工具箱与符号工具箱的应用,4.1曲线拟合与插值运算,1多项式的建立与表示方法在MATLAB中，n次多项式用一个长度为n+1的行向量表示，其元素为多项式的系数，按降幂排列，缺少的幂次项系数为0。例如，多项式在MATLAB中用向量p=1-12025116表示。,2多项式的运算,(1)多项式的加减运算多项式的加减运算就是其所对应的系数向量的加减运算。相加减的多项式必须表示成相同的次数，如果次数不同，应该把低次的多项式不足的高次项用0补足。,(2)多项式的乘除运算命令w=conv(u,v)表示多项式u和v相乘，例如在MATLAB中输入u=1234，v=102030，c=conv(u,v)返回c=1040100160170120conv指令可以嵌套使用，如conv(conv(a,b),c)。命令q,r=deconv(v,u)表示u整除v。向量q表示商，向量r表示余，即有v=conv(u,q)+r。,(3)多项式的导函数对多项式求导数的函数有k=polyder(p)，返回多项式p的导函数；k=polyder(a,b)，返回多项式a与b的乘积的导函数；q,d=polyder(b,a)，返回多项式b整除a的导函数，其分子多项式返回给q，分母多项式返回给d。,(4)多项式求值MATLAB中提供了两种求多项式值的函数。y=polyval(p,x)，代数多项式函数求值，若x为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。Y=polyvalm(p,x)，矩阵多项式求值，要求x为方阵。设A为方阵，p代表多项式x3-5x2+8，那么polyvalm(p,A)的含义是A*A*A-5*A*A+8*eye(size(A)而polyval(p,A)的含义是A.*A.*A-5*A.*A+8*ones(size(A),例4.1多项式P=x4-29x3+72x2-29x+1，以4阶pascal矩阵为自变量分别用polyval和polyvalm计算该多项式的值。在命令窗口输入如下命令：p=1-2972-291;X=pascal(4);A=polyval(p,X),B=polyvalm(p,X),返回A=161616161615-140-56316-140-2549-1208916-563-12089-43779B=0000000000000000,(5)多项式的根使用函数roots可以求出多项式等于0的根，根用列向量表示，其调用格式为r=roots(p)若已知多项式等于0的根，函数poly可以求出相应多项式，调用格式为p=poly(r),例4.2求多项式x4+8x3-10的根。命令如下：A=1,8,0,0,-10;x=roots(A)返回x=-8.01941.0344-0.5075+0.9736i-0.5075-0.9736i再输入p=poly(x)返回p=1.00008.0000-0.0000-0.0000-10.0000,3曲线拟合,在MATLAB中用polyfit函数来求得最小二乘拟合多项式的系数，再用polyval函数按所得的多项式计算所给出的点上的函数近似值。polyfit函数的调用格式为P,S=polyfit(X,Y,m)函数根据采样点X和采样点函数值Y，产生一个m次多项式P及其在采样点的误差向量S。其中X,Y是两个等长的向量，P是一个长度为m+1的向量，P的元素为多项式系数。,例4.3用一个6次多项式在区间0,2内逼近函数,MATLAB程序如下：x=linspace(0,2*pi,50);y=sin(x);P=polyfit(x,y,6)%得到6次多项式的系数和误差程序运行结果如下：P=0.0000-0.00560.0874-0.39460.26850.87970.0102,绘出sinx和多项式P(x)在给定区间的函数曲线，如图4.1所示。,图4.1用6次多项式对正弦函数进行拟合,4数据插值,(1)一维数据插值在MATLAB中，实现一维数据插值的函数是interp1，被插值函数是一个单变量函数，其调用格式为yi=interp1(x,Y,xi,method)函数根据x,y的值，计算函数在xi处的值。x,y是两个等长的已知向量，分别描述采样点和样本值，xi是一个向量或标量，描述欲插值的点，yi是一个与xi等长的插值结果。method是插值方法，允许的取值有linear、nearest、cubic、spline，分别表示线性插值、最近点插值、3次多项式插值、3次样条插值。注意：xi的取值范围不能超出x的给定范围，否则，会给出“NaN”错误。,例4.4向量t和p表示从19001990年的每隔10年的美国人口普查数据：t=1900:10:1990;p=75.99591.972105.711123.203131.669.150.697179.323203.212226.505249.633;根据人口普查数据估计1975年的人口，并利用插值估计从19002000年每年的人口数。,首先在命令窗口输入插值命令interp1(t,p,1975)，估计1975年的人口，返回结果如下：ans=214.8585再利用插值估计19002000年每年的人口数，并利用画图命令得出曲线分布。在MATLAB命令窗口输入如下命令语句：x=1900:1:2000;y=interp1(t,p,x,spline);%样条插值plot(t,p,:o,x,y,-r)%带圆圈标记的虚线(:o)为普查数据，红实线%(-r)为插值数据所得曲线如图4.2所示。,图4.2用一维数据插值得到的美国100年人口分布图,(2)二维数据插值在MATLAB中，提供了解决二维插值问题的函数interp2，其调用格式为Z1=interp2(X,Y,Z,X1,Y1,method)其中，X,Y是两个向量，分别描述两个参数的采样点；Z是与参数采样点对应的函数值；X1,Y1是两个向量或标量，描述欲插值的点；Z1是根据相应的插值方法得到的插值结果。method的取值与一维插值函数相同。X,Y,Z也可以是矩阵形式。同样，X1,Y1的取值范围不能超出X,Y的给定范围，否则，会给出“NaN”错误。,例4.5利用插值运算对峰值函数peaks插入更多的栅格。MATLAB程序如下：X,Y=meshgrid(-3:.25:3);Z=peaks(X,Y);XI,YI=meshgrid(-3:.125:3);ZI=interp2(X,Y,Z,XI,YI);mesh(X,Y,Z),hold,mesh(XI,YI,ZI+15)holdoffaxis(-33-33-520)程序运行结果如图4.3所示。,图4.3用二维数据插值得到的peaks函数曲线图,例4.6给定雇员数据如下：years=1950:10:1990;%工作年份service=10:10:30;%服役时间，即在职时间wage=150.697199.592187.625179.323195.072250.287203.212179.092322.767226.505153.706426.730249.633120.281598.243;%工资利用二维数据插值计算一个雇员在工作15年后在1975年所获得的工资。,MATLAB程序如下：w=interp2(service,years,wage,15,1975)程序运行结果如下：w=190.6287,4.2.1数值微分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时，按行计算差分。,4.2数值微积分,例4.7设x由0,2间均匀分布的10个点组成，求sinx的13阶差分。MATLAB程序如下：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),程序运行结果如下：X=00.69811.39632.09442.79253.49074.18884.88695.58516.2832Y=00.64280.98480.86600.3420-0.3420-0.8660-0.9848-0.6428-0.0000DY=0.64280.3420-0.1188-0.5240-0.6840-0.5240-0.11880.34200.6428D2Y=-0.3008-0.4608-0.4052-0.16000.16000.40520.46080.3008D3Y=-0.16000.05560.24520.32010.24520.0556-0.1600,4.2.2数值积分,quad(filename,a,b,tol,trace)quadl(filename,a,b,tol,trace)其中filename是被积函数名。a和b分别是定积分的下限和上限。tol用来控制积分精度，默认时为10-6。trace控制是否展现积分过程，若取非0则展现积分过程，若取0则不展现，默认时取trace=0。,例4.9用两种方法求在MATLAB命令窗口输入F=inline(1./(x.3-2*x-5);Q=quad(F,0,2)或Q=quadl(F,0,2),返回Q=-0.4605或者采用函数句柄Q=quad(myfun,0,2)也可以返回Q=-0.4605,这里myfun.m为一个M文件：functiony=myfun(x)y=1./(x.3-2*x-5);,二重定积分,I=dblquad(f,a,b,c,d,tol,trace)该函数求f在a,bc,d区域上的二重定积分。参数tol,trace的用法与函数quad完全相同。,例4.11计算二重定积分首先建立一个函数文件fxy.m，functionf=fxy(x,y)globalki;ki=ki+1;%ki用于统计被积函数的调用次数f=exp(-x.2/2).*sin(x.2+y);,再调用dblquad函数求解。globalki;ki=0;I=dblquad(fxy,-2,2,-1,1)ki程序运行结果如下：I=1.57449318974494ki=1038,4.3线性方程组求解,对于线性方程组Ax=b，可以利用左除运算符“”求解：x=Ab即x=inv(A)*b如果矩阵A是奇异的或接近奇异的，则MATLAB会给出警告信息。,MATLAB提供的lu函数用于对矩阵进行LU分解，其调用格式如下。L,U=lu(X)，产生一个上三角阵U和一个变换形式的下三角阵L(行交换)，使之满足X=LU。注意，这里的矩阵X必须是方阵。L,U,P=lu(X)，产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PX=LU。当然矩阵X同样必须是方阵。实现LU分解后，线性方程组Ax=b的解可表示为x=U(Lb)或x=U(LPb)，这样可以大大提高运算速度。,MATLAB的函数qr可用于对矩阵进行QR分解，其调用格式如下。Q,R=qr(X)，产生一个正交矩阵Q和一个上三角矩阵R，使之满足X=QR。Q,R,E=qr(X)，产生一个正交矩阵Q、一个上三角矩阵R以及一个置换矩阵E，使之满足XE=QR。实现QR分解后，线性方程组Ax=b的解可表示为x=R(Qb)或x=E(R(Qb)。,MATLAB函数chol(X)用于对矩阵X进行Cholesky分解，其调用格式如下。R=chol(X)，产生一个上三角阵R，使RR=X。若X为非对称正定的，则输出一个出错信息。R,p=chol(X)，这个命令格式将不输出出错信息。当X为对称正定的，则p=0，R与上述格式得到的结果相同；否则p为一个正整数。如果X为满秩矩阵，则R为一个阶数为q=p-1的上三角阵，且满足RR=X(1:q,1:q)。实现Cholesky分解后，线性方程组Ax=b变成RRx=b，所以x=R(Rb)。,例4.12用直接法求解下列线性方程组。,MATLAB程序如下：A=1/21/31/4;1/31/41/5;1/41/51/6B=0.950.670.52X=AB或X=inv(A)*B程序运行结果如下：X=1.20000.60000.6000,例4.13分别用LU分解、QR分解和Cholesky分解求解例4.12中的线性方程组。,MATLAB程序如下：A=1/21/31/4;1/31/41/5;1/41/51/6B=0.950.670.52L,U=lu(A);X=U(LB)%LU分解求解Q,R=qr(A);X=R(QB)%QR分解求解R=chol(A)X=R(RB)%Cholesky分解求解,所得结果均为X=1.20000.60000.6000,Jacobi迭代法的MATLAB函数文件Jacobi.m如下：functiony,n=jacobi(A,b,x0,eps)ifnargin=3eps=1.0e-6;elseifnargin=epsx0=y;y=B*x0+f;n=n+1;end,例4.14用Jacobi迭代法求解线性方程组。设迭代初值为0，迭代精度为10-6。,在命令中调用函数文件Jacobi.m，命令如下：A=9,-1,1;-1,10,-2;-2,1,10;b=10,7,6;x,n=jacobi(A,b,0,0,0,1.0e-6),MATLAB返回x=1.13650.95990.7313n=10输入x=inv(A)*b，所得结果相同。,4.4常微分方程的数值求解,MATLAB提供了求常微分方程数值解的函数，一般调用格式如下：t,y=ode23(fname,tspan,y0),二阶、三阶龙格-库塔法t,y=ode45(fname,tspan,y0)，四阶、五阶龙格-库塔法其中，fname是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量；tspan形式为t0,tf；表示求解区间；y0是初始状态列向量；t和y分别给出时间向量和相应的状态向量。,例4.17考虑著名的Rossler微分方程组,选定a=b=0.2，c=5.7，且x(0)=y(0)=z(0)=0,求解该微分方程。,引入新状态变量,其矩阵形式为,编写M函数如下：functiondx=rossler(t,x,flag,a,b,c)%虽然不显含时间，但还应该写出占位dx=-x(2)-x(3);x(1)+a*x(2);b+(x(1)-c)*x(3);%对应方程,保存为rossler.m文件，在MATLB命令窗口输入如下语句：a=0.2;b=0.2;c=5.7;%从函数外部定义参数变量x0=000;%微分方程的初值t,y=ode45(rossler,0,100,x0,a,b,c)%求解微分方程plot(t,y)%绘制各个状态变量的时间响应figureplot3(y(:,1),y(:,2),y(:,3),上面的命令直接得出该微分方程在内的数值解，并绘制出各个状态变量和时间之间的关系曲线和相空间曲线，如图4.5(a)和图4.5(b)所示。,图4.5Rossler方程的数值解表示,4.5MATLAB符号计算,(1)sym函数sym函数用来建立单个符号量，其调用格式为符号变量名=sym(符号字符串)符号字符串可以是常量、变量、函数或表达式。,(2)syms函数syms函数的一般调用格式为symsvar1var2varn函数定义符号变量var1,var2,varn等。用这种格式定义符号变量时不要在变量名上加字符分界符()，变量间用空格而不要用逗号分隔。symsxbeta等同于x=sym(x);beta=sym(beta);,建立符号表达式有以下两种方法。(1)用sym函数建立符号表达式命令语句y1=sym(1/sqrt(2*x)返回y1=1/sqrt(2*x)命令语句M=sym(a,b;c,d)返回M=a,bc,d,(2)使用已经定义的符号变量组成符号表达式命令语句symsxyv=3*x2-5*y+2*x*y+6返回v=3*x2-5*y+2*x*y+6,符号表达式的加、减、乘、除运算可分别由函数symadd，symsub，symmul和symdiv来实现，幂运算可以由sympow来实现。此外，和数值运算一样，也可以用+、-、*、/、运算实现符号运算。,n,d=numden(s)该函数提取符号表达式s的分子和分母，分别将它们存放在n与d中。例如：命令语句n,d=numden(sym(4/5)返回n=4andd=5,factor(S)，对S分解因式，S是符号表达式或符号矩阵。expand(S)，对S进行展开，S是符号表达式或符号矩阵。collect(S)，对S合并同类项，S是符号表达式或符号矩阵。collect(S,v)，对S按变量v合并同类项，S是符号表达式或符号矩阵。,例如：命令语句f=factor(123)返回f=341命令语句expand(x-2)*(x-4)返回x2-6*x+8命令语句expand(cos(x+y)返回cos(x)*cos(y)-sin(x)*sin(y)命令语句expand(exp(a+b)2)返回exp(a2)*exp(a*b)2*exp(b2),simplify(S)，应用函数规则对S进行化简。simple(S)，调用MATLAB的其他函数对表达式进行综合化简，并显示化简过程。例如：命令语句simplify(sin(x)2+cos(x)2)返回1命令语句simplify(exp(c*log(sqrt(a+b)返回(a+b)(1/2c)命令语句S=(x2+5*x+6)/(x+2),sqrt(16);R=simplify(S)返回R=x+3,4,利用函数sym可以将数值表达式变换成它的符号表达式。例如：命令语句sym(1.5)返回ans=3/2命令语句sym(3.14)返回ans=157/50,函数numeric或eval可以将符号表达式变换成数值表达式。N=numeric(S)等同于N=double(sym(S)，例如：命令语句A=magic(4);A(:,:,2)=A;d1,d2,d3=eval(size(A)返回d1=4d2=4d3=2,transpose(S)，返回S矩阵的转置矩阵。determ(S)，返回S矩阵的行列式值。colspace(S)，返回S矩阵列空间的基。,ezplot(f)，绘制函数f(x)在默认区间-2x2内的图形；ezplot(f,min,max)，绘制函数f(x)在指定区间min,max内的图形。该函数打开标签为FigureNo.1的图形窗口，并显示图形。如果已经存在图形窗口，该函数在标签数最大的窗口中显示图形；ezplot(f,xminxmax,fign)，在指定的窗口fign中绘制函数f(x)在指定区间min,max内的图形。,对于隐函数，ezplot函数的调用格式有：ezplot(f)，绘制函数f(x)在区间-2x2内的图形；ezplot(f,xmin,xmax,ymin,ymax)，绘制函数f(x,y)在指定区间xminxxmax，yminyymax的图形；ezplot(f,min,max)，绘制函数f(x,y)在指定区间minxmax，minymax的图形。,对于参数方程，ezplot函数的调用格式有：ezplot(x,y)，绘制参数方程x=x(t)，y=y(t)在默认区间0t2内的图形；ezplot(x,y,tmin,tmax)，绘制参数方程x=x(t)，y=y(t)在指定区间tinttmax的图形。,ezplot3函数用于绘制三维参数曲线。该函数的调用格式有：ezplot3(x,y,z)，在默认区间0t2内绘制参数方程x=x(t)，y=y(t)，z=z(t)的图形；ezplot3(x,y,z,tmin,tmax)，在区间tminttmax内绘制参数方程x=x(t)，y=y(t)，z=z(t)的图形；ezplot3(.,animate)，生成空间曲线的动态轨迹。,ezmesh，ezsurf函数分别用于绘制三维网格曲面图和三维表面图。这两个函数的用法相同，ezmesh函数的调用格式有：ezmesh(f)，绘制函数f(x,y)在默认区间-2x2，-2y2内的图形；ezmesh(f,domain)，在指定区域domain绘制函数f(x,y)的图形；ezmesh(x,y,z)，绘制参数方程x=x(s，t)，y=y(s，t)，z=y(s，t)在默认区间-2s2，-2t2内的图形,ezmesh(x,y,z,smin,smax,tmin,tmax)或ezmesh(x,y,z,min,max)，在指定区域绘制三维参数方程的图形。ezmeshc，ezsurfc两个函数用于在绘制三维曲面的同时绘制等值线。ezmeshc(f)，绘制二元函数f(x,y)在默认区间-2x2，-2y2的图形。ezmeshc(f,domain)，绘制函数f(x,y)在指定区域的图形，绘图区域由domain指定，其中domain为41数组或者21数组，如xmin,xmax,ymin,ymax表示minxmax，minymax，min,max表示minxmax,minymax。,ezmeshc(x,y,z)，绘制参数方程x=x(s,t)，y=y(s,t)和z=z(s,t)在默认区间-2s2，-2t2的图形。ezmeshc(x,y,z,smin,smax,tmin,tmax)，ezmeshc(x,y,z,min,max)，绘制参数方程在指定区域的图形。ezmeshc(.,n)，指定绘图的网格数，默认值为60。ezmeshc(.,circ)，在以指定区域中心为中心的圆盘上绘制图形,在MATLAB中，用于绘制符号函数等值线的函数有ezcontour和ezcontourf，这两个函数分别用于绘制等值线和带有区域填充的等值线。其调用格式有：ezcontour(f)，绘制符号二元函数f(x,y)在默认区间-2x2，-2y2内的等值线图；ezcontour(f,domain)，绘制符号二元函数f(x,y)在指定区域的等值线图；ezcontour(.,n)，绘制等值线图，并指定等值线的数目。例如命令ezplot(x2-y4)绘制的图形曲线如图4.7所示。例如命令ezmeshc(y/(1+x2+y2),-5,5,-2*pi,2*pi)绘制的图形如图4.8所示。,limit(f,x,a)，计算当变量x趋近于常数a时，f(x)函数的极限值。limit(f,a)，求符号函数f(x)的极限值，符号函数f(x)的变量为函数findsym(f)确定的默认自变量，即变量x趋近于a。limit(f)，系统默认变量趋近于0，即a=0的极限。limit(f,x,a,right)，变量x从右边趋近于a时符号函数f(x)的极限值。limit(f,x,a,left)，变量x从左边趋近于a时符号函数的极限值。,diff(f,x,n)diff函数求函数f对变量x的n阶导数。参数x的用法同求极限函数limit，可以默认，默认值与limit相同，n的默认值是1。,int(f,x)，求函数f对变量x的不定积分。参数x可以默认，默认原则与diff函数相同。int(f,v)，以v为自变量，对被积函数符号表达式f求不定积分。int(f,v,a,b)，求被积函数f在区间a,b上的定积分。a和b可以是两个具体的数，也可以是一个符号表达式。,例4.20求椭球的体积。,MAT