实验数学五：MATLAB在数值计算中的应用

1MATLAB在数值计算中的应用5.1实验目的在工程技术中，大量的实际问题都需要进行近似处理，从而产生不同问题的数值计算方法。而MATLAB具有强大的数值运算功能，本实验的目的是学会用MATLAB软件进行一些数值运算，包括代数方程求根、插值问题和曲线拟合问题等。5.2实验内容一、代数方程求根代数方程求根有各种近似处理方法，下面给出MATLAB两种常用的调用格式：最小二乘法格式：fsolve(‘f’,x0):求方程f=0在估计值x0附近的近似解。2例100x-x求方程x-e在附近的根。解输入命令：f=inline('x-exp(-x)');x1=fsolve(f,0)x1=0.5671例225x-x求方程sinx-e=0的根。先画图观察根的个数及大概位置。输入命令：fplot('[5*x^2*sin(x)-exp(-x),0]',[0,10])结果见图5.1注意，[5*x^2*sin(x)-exp(-x),0]中的[…,0]是作y=0直线，即x轴。3方程在[0,10]区间从图中可看出有4个解，分别在0，3，6，9附近，所以用命令：f=inline('5*x.^2.*sin(x)-exp(-x)');fsolve(f,[0,3,6,9])ans=0.50183.14076.28329.424842、零点法格式：fzero(‘f’,x0):求函数f在x0附近的零点。例32求方程x-4x-5=0的根。先画图观察根的个数及大概位置。输入命令：fplot('[x^2-4*x-5,0]',[-10,10])结果见图5.2fzero(‘f’,[x1,x2]):求函数f在区间[x1,x2]上唯一零点。5从图中可看出方程在[-2,0]及[4，6]区间上各有一根，再输入命令：x1=fzero('x^2-4*x-5',[-2,0])x1=-1x2=fzero('x^2-4*x-5',[4,6])x2=563、代数方程的符号解格式：solve('f',):求代数方程f=0的根；solve('eqn1','eqn2',...,'eqnN'):求n个代数方程的根；例40bxc2求方程ax的解解输入命令：solve('a*x^2+b*x+c')ans=[1/2/a*(-b+(b^2-4*a*c)^(1/2))][1/2/a*(-b-(b^2-4*a*c)^(1/2))]7例5115xyxy=1求方程组解输入命令：[x,y]=solve('x*y=1','x-11*y=5')x=[5/2+1/2*69^(1/2)][5/2-1/2*69^(1/2)]y=[-5/22+1/22*69^(1/2)][-5/22-1/22*69^(1/2)]如果化成数值解，用命令vpa如上例：x=vpa(x,2)x=[6.7][-1.7]8y=vpa(y,2)y=[.14][-.60]二、曲线拟合已知离散点上的数据集1122[(,)(,)(,)],nnxyxyxy求得一解析函数y=f(x)使y=f(x)在原离散点ix接近给定iy曲线拟合是最小二乘法曲线拟合，拟合结果可使误差的上尽可能的值，这一过程叫曲线拟合。最常用的平方和最小，即找出使21()niiifxy最小的f(x).9格式：p=polyfit(x,y,n).说明：求出已知数据x,y的n次拟合多项式f(x)的系数p,x必须是单调的。例6已知某函数的离散值如表5.1xi0.51.01.52.02.53.0yi1.752.453.814.807.008.65求二次拟合多项式.先画函数离散点的图形输入命令：x=[0.51.01.52.02.53.0];y=[1.752.453.814.807.008.60];scatter(x,y,5)结果见图5.310由图可看出可用二次多项式拟合。再输入命令：p=polyfit(x,y,2)p=0.56140.82871.1560即二次拟合多项式为2()0.56140.82871.1560fxxx11画出离散点及拟合曲线：输入命令：x1=0.5:0.05:3.0;y1=polyval(p,x1);plot(x,y,'*r',x1,y1,'-b')结果见图5.412三、一维插值已知离散点上的数据集1122[(,)(,)(,)],nnxyxyxy求得一解析函数连接自变量相邻的两个点，并求得两点间的数值，这一过程叫插值。MATLAB在一维插值函数interp1中，提供了四种插值方法选择：线性插值、三次样条插值、立方插值和最近邻点插值。interp1的本格式为：yi=interp1(x，y，xi，'method')其中x,y分别表示数据点的横、纵坐标向量，x必须单调，xi为需要插值的横坐标数据（或数组），xi不能超出x的范围，而method为可选参数，有四种选择：‘nearest’：最邻近插值‘linear’：线性插值；13‘spline’：三次样条插值；‘cubic’：立方插值。缺省时：分段线性插值。例7在1-12的11小时内，每隔1小时测量一次温度，测得的温度依次为：5，8，9，15，25，29，31，30，22，25，27，24。试估计在3.2,6.5,7.1,11.7小时的温度值。解输入命令：hours=1:12;temps=[589152529313022252724];t=interp1(hours,temps,[3.26.57.111.7])%线性插值t=10.200030.000030.900024.900014T=interp1(hours,temps,[3.26.57.111.7],'spline')%三次样条插值T=9.673430.042731.175525.3820比较发现，两种结果有差异，这是因为插值是一个估计或猜测的过程。两种插值的画图如下；输入命令：t0=1:0.1:12;T0=interp1(hours,temps,t0,'spline');plot(hours,temps,'+',t0,T0,hours,temps,'r:')xlabel('时间');ylabel('温度')15gtext('线性插值')gtext('三次样条插值')结果见图5.516四、二维插值对二维插值问题，MATLAB分别给出了针对插值基点为网格节点的插值函数及针对插值基点为散乱节点的插值函数调用格式。1、用MATLAB作网格节点数据的插值;1,2,,mjnijij已知mn个节点：(x,y,z),(i=1,2,)12;nxxx且***(,)((,)).ijxyxyz求点处的插值12,nyyy对上述问题，MATLAB提供了二维插值函数interp2,其基本格式为：17z=interp2(x0,y0,z0,x,y,’method’)其中x0,y0是自变量。X0,y0的分量值必须是单调递增的。X0和y0分别是m维和n维向量，分别表示已知数据点的横、纵坐标向量，z0是m*n维矩阵，标明相应于所给数据网格点的函数值。向量x,y是待求函数值所给定网格点的的横、纵坐标向量，x,y的值分别不能超出x0,y0的范围。而method为可选参数，有四种选择：‘nearest’最邻近插值‘linear’线性插值‘spline’三次样条插值‘cubic’三次插值缺省时,是线性插值18例8：测得平板表面3×5网格点处的温度分别为：828180828479636165818484828586试求在平板表面坐标为（1.5,1.5）,（2,1.6）,（2.5,2）（3.5,4.5）处的温度，并作平板表面的温度分布曲面z=f(x,y)的图形，（1）先在三维坐标画出原始数据，画出粗糙的温度分布曲图.输入以下命令：x0=1:5;y0=1:3;temps=[8281808284;7963616581;8484828586];19mesh(x0,y0,temps)结果见图5.6分别用线性性插值和三次样条插值求已知点的温度。输入命令：20t=interp2(x0,y0,temps,[1.522.53.5],[1.51.624.5],'liner')t=76.250070.200062.0000NaNT=interp2(x0,y0,temps,[1.522.53.5],[1.51.624.5],'spline')T=71.453165.520060.9688188.8906（2）以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值画出线性和三次样条插值的温度分布曲面图.输入以下命令得温度的线性插值曲面图:x=1:0.2:5;y=1:0.2:3;z=interp2(x0,y0,temps,x',y,'linear');mesh(x,y,z)21xlabel('x轴');ylabel('y轴');zlabel('z轴')title('线性插值曲面图')结果见图5.722再输入以下命令得温度的三次样条插值曲面图:z=interp2(x0,y0,temps,x',y,'spline');mesh(x,y,z)xlabel('x轴');ylabel('y轴');zlabel('z轴')title('三次样条插值曲面图')结果见图5.723例9山区地貌：在某山区测得一些地点的高程如下表5.2。平面区域为1200=x=4000,1200=y=3600)试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。Yx04008001200160020002400280032003600400044004800520056000400800120016002000240028003200360040004400480037047055060067069067062058045040030010015025051062073080085087085078072065050020030035032065076088097010201050102083090070030050055048035074088010801130125012801230104090050070078075065055083098011801320145014201400130070090085084038078075088010601230139015001500140090011001060950870900930950910109012701500120011001350145012001150101088010001050110095011901370150012001100155016001550138010709001050115012001430143014601500155016001550160016001600155015001500155015501420143014501480150015501510143013001200980850750550500138014101430145014701320128012001080940780620460370350137013901410143014401140111010509508206905403803002101350137013901400141096094088080069057043029021015024此例将对最近邻点插值、线性插值方法和双三次插值方法的插值效果进行比较。输入命令：x0=0:400:5600;y0=0:400:4800;z0=[370470550600670690670620580450400300100150250;...510620730800850870850780720650500200300350320;...650760880970102010501020830900700300500550480350;...