matlab数值计算与符号计算

第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]表示。43212025116xxxx2．多项式的运算•(1)多项式的加减运算•多项式的加减运算就是其所对应的系数向量的加减运算。相加减的多项式必须表示成相同的次数，如果次数不同，应该把低次的多项式不足的高次项用0补足。•(2)多项式的乘除运算•命令w=conv(u,v)表示多项式u和v相乘，例如在MATLAB中输入•u=[1234]，v=[102030]，c=conv(u,v)•返回•c=•1040100160170120•conv指令可以嵌套使用，如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=•16161616•1615-140-563•16-140-2549-12089•16-563-12089-43779•B=•0000•0000•0000•0000•(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.0194•1.0344•-0.5075+0.9736i•-0.5075-0.9736i•再输入•p=poly(x)•返回•p=•1.00008.0000-0.0000-0.0000-10.00003．曲线拟合•在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表示从1900～1990年的每隔10年的美国人口普查数据：•t=1900:10:1990;•p=[75.99591.972105.711123.203131.669...•150.697179.323203.212226.505249.633];•根据人口普查数据估计1975年的人口，并利用插值估计从1900～2000年每年的人口数。•首先在命令窗口输入插值命令interp1(t,p,1975)，估计1975年的人口，返回结果如下：•ans=•214.8585•再利用插值估计1900～2000年每年的人口数，并利用画图命令得出曲线分布。在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)•holdoff•axis([-33-33-520])•程序运行结果如图4.3所示。•图4.3用二维数据插值得到的peaks函数曲线图•例4.6给定雇员数据如下：•years=1950:10:1990;%工作年份•service=10:10:30;%服役时间，即在职时间•wage=[150.697199.592187.625•179.323195.072250.287•203.212179.092322.767•226.505153.706426.730•249.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的1～3阶差分。•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.2832•Y=•00.64280.98480.86600.3420-0.3420-0.8660-0.9848-0.6428-0.0000•DY=•0.64280.3420-0.1188-0.5240-0.6840-0.5240-0.11880.34200.6428•D2Y=•-0.3008-0.4608-0.4052-0.16000.16000.40520.46080.3008•D3Y=•-0.16000.05560.24520.32010.24520.0556-0.16004.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)2301d25Ixxx•返回•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,b]×[c,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);212/2212esin()ddxIxyxy•再调用dblquad函数求解。•globalki;ki=0;•I=dblquad('fxy',-2,2,-1,1)•ki•程序运行结果如下：•I=•1.57449318974494•ki=•10384.3线性方程组求解•对于线性方程组Ax=b，可以利用左除运算符“\”求解：•x=A\b•即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\(L\