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

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,MATLAB在数值计算中的应用,5.1实验目的,在工程技术中，大量的实际问题都需要进行近似处理，从而产生不同问题的数值计算方法。而MATLAB具有强大的数值运算功能，本实验的目的是学会用MATLAB软件进行一些数值运算，包括代数方程求根、插值问题和曲线拟合问题等。,5.2实验内容,一、代数方程求根,代数方程求根有各种近似处理方法，下面给出MATLAB两种常用的调用格式：,最小二乘法,格式：fsolve(f,x0):求方程f=0在估计值x0 附近的近似解。,1,例1,解,输入命令 ：, f=inline(x-exp(-x);, x1=fsolve(f,0),x1 =,0.5671,例2,先画图观察根的个数及大概位置。,输入命令 ：, fplot(5*x2*sin(x)-exp(-x),0,0,10),结果见图5.1,注意，5*x2*sin(x)-exp(-x),0中的 ,0是作y=0直线，即x轴。,2,方程在0,10区间从图中可看出有4个解，分别在0，3，6，9附近，,所以用命令：, f=inline(5*x.2.*sin(x)-exp(-x);, fsolve(f,0,3,6,9),ans =,0.5018 3.1407 6.2832 9.4248,3,2、零点法,格式：fzero(f,x0): 求函数f在x0 附近的零点。,例3,先画图观察根的个数及大概位置。,输入命令 ：, fplot(x2-4*x-5,0,-10,10),结果见图5.2,fzero(f,x1,x2): 求函数f在区间x1,x2 上唯一零点。,4,从图中可看出方程在-2,0及4，6区间上各有一根，,再输入命令 ：, x1=fzero(x2-4*x-5,-2,0),x1 =,-1, x2=fzero(x2-4*x-5,4,6),x2 =,5,5,3、代数方程的符号解,格式：,solve(f,):求代数方程f=0的根；,solve(eqn1,eqn2,.,eqnN):求n个代数方程的根；,例4,解,输入命令 ：, solve(a*x2+b*x+c),ans =, 1/2/a*(-b+(b2-4*a*c)(1/2), 1/2/a*(-b-(b2-4*a*c)(1/2),6,例5,解,输入命令 ：, 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,7, y=vpa(y,2),y =, .14, -.60,二、曲线拟合,已知离散点上的数据集,求得一解析函数y=f(x)使y=f(x)在原离散点,接近给定,曲线拟合是最小二乘法曲线拟合，拟合结果可使误差的,上尽可能,的值，这一过程叫曲线拟合。最常用的,平方和最小，即找出使,最小的f(x).,8,格式：p=polyfit(x,y,n).,说明：求出已知数据x,y 的n次拟合多项式f(x)的系数p,x 必须是单调的。,例6 已知某函数的离散值如表5.1,xi,0.5,1.0,1.5,2.0,2.5,3.0,yi,1.75,2.45,3.81,4.80,7.00,8.65,求二次拟合多项式.,先画函数离散点的图形,输入命令 ：, x=0.5 1.0 1.5 2.0 2.5 3.0;, y=1.75 2.45 3.81 4.80 7.00 8.60;, scatter(x,y,5),结果见图5.3,9,由图可看出可用二次多项式拟合。,再,输入命令 ：, p=polyfit(x,y,2),p =,0.5614 0.8287 1.1560,即二次拟合多项式为,10,画出离散点及拟合曲线：,输入命令 ：, x1=0.5:0.05:3.0;, y1=polyval(p,x1);, plot(x,y,*r,x1,y1,-b),结果见图5.4,11,三、一维插值,已知离散点上的数据集,求得一解析函数连接自变量相邻的两个点，并求得两点,间的数值，这一过程叫插值。,MATLAB在一维插值函数interp1中，提供了四种插值方法选择：线性插值、三次样条插值、立方插值和最近邻点插值。interp1的本格式为：,yi=interp1(x，y，xi，method),其中x,y分别表示数据点的横、纵坐标向量，x 必须单调，xi为需要插值的横坐标数据（或数组），xi不能超出x的范围，而method为可选参数，有四种选择：,nearest ：最邻近插值,linear ： 线性插值；,12,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=5 8 9 15 25 29 31 30 22 25 27 24;, t=interp1(hours,temps,3.2 6.5 7.1 11.7) %线性插值,t =,10.2000 30.0000 30.9000 24.9000,13, T=interp1(hours,temps,3.2 6.5 7.1 11.7,spline) %三次样条插值,T =,9.6734 30.0427 31.1755 25.3820,比较发现，两种结果有差异，这是因为插值是一个估计或猜测的过程。,两种插值的画图如下；,输入命令 ：, t0=1:0.1:12;, T0=interp1(hours,temps,t0,spline);, plot(hours,temps,+,t0,T0,hours,temps,r:), xlabel(时间);, ylabel(温度),14, gtext(线性插值), gtext(三次样条插值),结果见图5.5,15,四、二维插值,对二维插值问题，MATLAB分别给出了针对插值基点为网格节点的插值函数及针对插值基点为散乱节点的插值函数调用格式。,1、,用MATLAB作网格节点数据的插值,对上述问题，MATLAB提供了二维插值函数interp2,其基本格式为：,16,z=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 三次插值,缺省时, 是线性插值,17,例8：测得平板表面35网格点处的温度分别为：,82 81 80 82 84 79 63 61 65 81 84 84 82 85 86,试求在平板表面坐标为（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=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;,18,mesh(x0,y0,temps),结果见图5.6,分别用线性性插值和三次样条插值求已知点的温度。,输入命令 ：,19, t=interp2(x0,y0,temps,1.5 2 2.5 3.5,1.5 1.6 2 4.5,liner),t =,76.2500 70.2000 62.0000 NaN, T=interp2(x0,y0,temps,1.5 2 2.5 3.5,1.5 1.6 2 4.5,spline),T =,71.4531 65.5200 60.9688 188.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),20, xlabel(x轴);, ylabel(y轴);, zlabel(z轴), title(线性插值曲面图),结果见图5.7,21,再输入以下命令得温度的三次样条插值曲面图:, z=interp2(x0,y0,temps,x,y,spline);, mesh(x,y,z), xlabel(x轴);, ylabel(y轴);, zlabel(z轴), title(三次样条插值曲面图),结果见图5.7,22,例 9 山区地貌：,在某山区测得一些地点的高程如下表5.2。平面区域为 1200=x=4000,1200=y x0=0:400:5600;, y0=0:400:4800;,z0=370 470 550 600 670 690 670 620 580 450 400 300 100 150 250;.,510 620 730 800 850 870 850 780 720 650 500 200 300 350 320;.,650 760 880 970 1020 1050 1020 830 900 700 300 500 550 480 350;.,740 880 1080 1130 1250 1280 1230 1040 900 500 700 780 750 650 550;.,830 980 1180 1320 1450 1420 1400 1300 700 900 850 840 380 780 750;.,880 1060 1230 1390 1500 1500 1400 900 1100 1060 950 870 900 930 950;.,910 1090 1270 1500 1200 1100 1350 1450 1200 1150 1010 880 1000 1050 1100;.,950 1190 1370 1500 1200 1100 1550 1600 1550 1380 1070 900 1050 1150 1200;.,1430 1430 1460 1500 1550 1600 1550 1600 1600 1600 1550 1500 1500 1550 1550;.,1420 1430 1450 1480 1500 1550 1510 1430 1300 1200 980 850 750 550 500;.,1380 1410 1430 1450 1470 1320 1280 1200 1080 940 780 620 460 370 350;.,1370 1390 1410 1430 1440 1140 1110 1050 950 820 690 540 380 300 210;.,1350 1370 1390 1400 1410 960 940 880 800 690 570 430 290 210 150;,24, meshz(x0,y0,z0), xlabel(x轴);, ylabel(y轴);, zlabel(z轴), title(原始数据山区表面图),结果见图5.8,25,每隔50加密网格，分别作最邻近插值、线性插值、三次样条插值，作出插值后的表面图：,（1）最近邻点插,输入命令 ：, x=0:50:5600;, y=0:50:4800;, z1=interp2(x0,y0,z0,x,y,nearest);, surfc(x,y,z1), xlabel(x轴);, ylabel(y轴);, zlabel(z轴), title(最近邻点插表面图),结果见图5.9,26,（2）线性插值,输入命令 ：, z2=interp2(x0,y0,z0,x,y,linear);, surfc(x,y,z2), xlabel(x轴);, ylabel(y轴);, zlabel(z轴), title(线性插值表面图),结果见图5.10,27,（3） 三次样条插值,输入命令 ：, z3=interp2(x0,y0,z0,x,y,spline);, surfc(x,y,z3), xlabel(x轴);, ylabel(y轴);, zlabel(z轴), title(三次样条插值表面图),结果见图5.11,28,（4）下面分别作出以上三种插值的等高线图,输入命令 ：, subplot(1,3,1),contour(x,y,z1,10,r);, subplot(1,3,2),contour(x,y,z2,10,r);, subplot(1,3,3),contour(x,y,z3,10,r);,结果见图5.12,从以上图形可以比较出三种插值法的效果。,29,2、,用MATLAB作散乱节点数据的插值,对上述问题，MATLAB提供了二维插值函数griddata,其基本格式为：,z =griddata（x0，y0，z0， x， y，method）,其中X0、y0、z0均是n维向量，分别表示已知数据点的横、纵坐标和竖坐标向量。向量x,y是待求函数值所给定网格点的的横、纵坐标向量。,而method为可选参数，有四种选择：,nearest 最邻近插值,linear 线性插值,cubic 三次插值,30,v4 Matlab提供的插值方法,缺省时, 线性插值,例10 在某海域测得一些点(x,y)处的水深z由下表5.3给出，船的吃水深度为5英尺，在矩形区域（75，200）*（-50，150）里的哪些地方船要避免进入。,X,129.0,140.0,103.5,88.0,185.5,195.0,105.5,157.5,107.5,77.0,81.0,162.0,162.0,117.5,Y,7.5,141.5,23.0,147.0,22.5,137.5,85.5,6.5,81.0,3.0,56.5,66.5,84.0,33.5,Z,4,8,6,8,6,8,8,9,9,8,8,9,4,9,31,(1)作出测量点的分布图,MATLAB命令：, clear;, x0=129 140 103.5 88 185.5 195 105.5 157.5 107.5 77 81 162 162 117.5;, y0=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5;, plot(x0,y0,+);,结果见图5.13,32,（2）作出三次插值海底地形图,输入命令 ：, x0=129 140 103.5 88 185.5 195 105.5 157.5 107.5 77 81 162 162 117.5;, y0=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5;, z0=-4 -8 -6 -8 -6 -8 -8 -9 -9 -8 -8 -9 -4 -9;, y=-70:0.5:150;, z=griddata(x0,y0,z0,x,y,cubic);, meshz(x,y,z);, xlabel(x轴);, ylabel(y轴);, zlabel(z轴), title(三次插值海底地形图), x=75:0.5:200;,33,结果见图5.14,34,（3）作出水深低于5英尺的海域范围,输入命令 ：, contour(x,y,z,-5,-5,-k);, grid, xlabel(x轴);, ylabel(y轴);, title(船不能进入区域),结果见图5.15,35,从图可以看出，船应避免进入危险区域,（-10，110）。,（110，190）,36,5.3 实验任务,3，5，6，8,37,