实验二 绘图及简单科学计算.doc

实验二 绘图及简单科学计算主要内容：1) 绘图简介2) 科学计算，曲线拟合与插值绘图简介MATLAB不但擅长于矩阵相关的数值运算，也适合用于各种科学可视化（Scientific visualization）。下面介绍MATLAB基本二维和三维的各项绘图命令，包含一维曲线及二维曲面的绘制、打印及保存。一、 二维绘图：l plot是绘制一维曲线的基本函数，但在使用此函数之前，我们需先定义曲线上每一点的x及y坐标。下例可画出一条正弦曲线： x=linspace(0, 2*pi, 100); % 100个点的x坐标 y=sin(x); % 对应的y坐标 plot(x,y); l 若要画出多条曲线，只需将坐标对依次放入plot函数即可： plot(x, sin(x), x, cos(x); l 若要改变颜色，在坐标对后面加上相关字符串即可： plot(x, sin(x), c, x, cos(x), g); l 若要同时改变颜色及线型（Line style），也是在坐标对后面加上相关字符串即可： plot(x, sin(x), co, x, cos(x), g*); plot绘图函数的参数字符颜色字符图线型态黄色.点k黑色o圆w白色xxb蓝色+g绿色*r红色-实线c亮青色:点线m锰紫色-.点虚线-虚线l 图形完成后，我们可用axis(xmin,xmax,ymin,ymax)函数来调整坐标轴的范围： axis(0, 6, -1.2, 1.2); l MATLAB也可对图形加上各种注解与处理： xlabel(Input Value); % x轴注解 ylabel(Function Value); % y轴注解 title(Two Trigonometric Functions); % 图形标题 legend(y = sin(x),y = cos(x); % 图形注解 grid on; % 显示格线 l 用subplot来同时画出数个小图形于同一个窗口之中： subplot(2,2,1); plot(x, sin(x); %把窗口分成2*2个子窗口,在第一个子窗口绘图 subplot(2,2,2); plot(x, cos(x); %在第二个子窗口绘图 subplot(2,2,3); plot(x, sinh(x); %在第三个子窗口绘图 %在第四个子窗口绘图l MATLAB还有其他各种二维绘图函数，以适合不同的应用，详见下表。 Bar长条图Errorbar图形加上误差范围Fplot较精确的函数图形Polar极坐标图Hist累计图Rose极坐标累计图Stairs阶梯图以下我们针对每个函数举例。l 当数据点数量不多时，条形图是很适合的表示方式： close all; % 关闭所有的图形窗口 x=1:10; y=rand(size(x); bar(x,y); 柱状图l 如果已知数据的误差量，就可用errorbar来表示。下例以单位标准差来做数据的误差量： x = linspace(0,2*pi,30); y = sin(x); e = std(y)*ones(size(x); errorbar(x,y,e) l 对于变化剧烈的函数，可用fplot来进行较精确的绘图，会对剧烈变化处进行较密集的取样，如下例： fplot(sin(1/x), 0.02 0.2); % 0.02 0.2是绘图范围 l 若要产生极坐标图形，可用polar： theta=linspace(0, 2*pi); r=cos(4*theta); polar(theta, r); l 对于大量的数据，我们可用hist来显示数据的分布情况和统计特性。下面几个命令可用来验证randn产生的高斯随随机数分布： x=randn(5000, 1); % 产生5000个 m=0，s=1 的高斯随机数 hist(x,20); % 20代表长条的个数 l rose和hist很接近，只不过是将数据大小视为角度，数据个数视为距离，并用极坐标绘制表示： x=randn(1000, 1); rose(x); l stairs可画出阶梯图： x=linspace(0,10,50); y= sin(x).*exp(-x/3); stairs(x,y); 1、 y=1 2 3 4;3 4 5 6;2 3 5 8; plot(y)2、 x=2 3 4 5; y=1 2 3 4;3 4 5 6;2 3 5 8; plot(x,y)3、 练习：还有其他可能情况，请同学们实验4、 hold on hold off 介绍5、 练习：利用已学知识绘制分段函数图像;6、 图轴控制：（1） 画出两个刻度不同的y 轴x=linspace(0,2*pi);y1=sin(x);y2=exp(-x);plotyy(x,y1,x,y2)（2） axis(xmin,xmax,ymin,ymax)(3) plot(x,y1,-,x,y2,:o) xlabel(x轴); ylabel(y轴)title(sin(x)与e-x比较)legend(sin(x),e-x) %多条曲线的说明 gtext(交点) grid on（3） 改变标刻点：x=0:0.1:4*pi;plot(x,sin(x)+sin(3*x)set(gca,ytick,-1,-0.3,0.1,1)set(gca,yticklabel极小,临界值,崩溃值,极大)grid on11、fplot,ezplot介绍fplot(sin(2*x)+cos(x),-10,10)ezplot(x2+1)12、介绍whos,clear all,clc,clf二、 三维图形 在科学可视化（Scientific visualization）中，三维空间的立体图是一个非常重要的技巧。下面介绍MATLAB基本三维空间的各项绘图命令。 l mesh和plot是空间立体绘图的基本命令，mesh可画出立体网状图，surf则可画出立体曲面图，两者产生的图形都会依高度而有不同颜色，surf和mesh的用法类似。l x=linspace(-2,2,25);y=linspace(-2,2,25);X,Y=meshgrid(x,y);%形成25阶的矩阵Z=X.*exp(-X.2-Y.2);mesh(X,Y,Z) %绘制立体网状图surf(X,Y,Z) %绘制立体曲面图l meshc同时画出网状图与等高线： x,y,z=peaks; meshc(x,y,z); axis(-inf inf -inf inf -inf inf); l surfc同时画出曲面图与等高线, 带等高线的贴面图： (1) x,y,z=peaks; surfc(x,y,z); axis(-inf inf -inf inf -inf inf); (2)clf;x=-8:0.5:8;y=x;X,Y=meshgrid(x,y); %用来产生三维绘图时的阵列（网格）Z=4*X.2-2*Y.2;surfc(X,Y,Z)l contour3画出曲面在三维空间中的等高线： contour3(peaks, 20); axis(-inf inf -inf inf -inf inf); l contour画出曲面等高线在XY平面的投影： c=contour(peaks,6)clabel(c) %给c阵所表示的等高线加注高度标识clabel(c,manual) %对鼠标击中的等高线加注高度标识注：由contourc所给出的等高线坐标阵c是一个行数为2的矩阵，它的排列方式是每一条由该矩阵的一个全行子阵表示，子阵的第一行第一列是等高线的高度值，子阵的第二行第一列是画该等高线的数据点数，子阵的其余列中的每一列都是一个数据点。l plot3可画出三维空间中的曲线, 也可同时画出两条三维空间中的曲线 t=linspace(0,20*pi, 501); plot3(t.*sin(t), t.*cos(t), t); t=linspace(0, 10*pi, 501); plot3(t.*sin(t), t.*cos(t), t, t.*sin(t), t.*cos(t), -t); (或者plot3(t.*sin(t), t.*cos(t), t,b-, t.*sin(t), t.*cos(t), -t,gd);)数组建立 字符串string=数学建模培训strin=e;e同一类型异质数组：其作用是可以储存不同的数据s=e,3,35 556 67建立内容显示celldisp(s)取用：s2科学计算多项式理论：多项式的表达：p=0.69552,0.436,0.668,1.35利用行向量表示降幂多项式的显示 poly2str(p,x)roots 求多项式的根 roots(p)poly：利用多项式的根，求出多项式的系数行向量，只需乘以原多项式最高次数项的系数 pp=poly(ans)，p=0.69552*ppfzero;找出函数值为零点，区间或起始点。fzero(x2+2*x,2)fzero(x2+2*x,-3,-1)Fminbnd;找出某区间的最小值fminbnd(x2-2*x-1+1/(x-0.8)2,-2,0.5)多项式的加减法：f(x)=1.35+0.668x+0.436x2+0.69552x3f(x)=2.3+0.453x+2.342x3p1=0.69552,0.436,0.668,1.35p2=2.342,0,0.45,2.3 p1+p2 %缺位补零多项式的乘法conv(p1,p2)多项式的除法;f(x)=1.35+0.668x+0.436x2+0.69552x3f(x)=1.23+0.645x+0.876x2p1=0.69552,0.436,0.668,1.35p2=0.876,0.645,1.23 a,b=deconv(p1,p2) %无需补零多项式的导数 polyder(p1)分式导数p1=0.69552,0.436,0.668,1.35p2=0.876,0.645,1.23polyder(p1,p2) a,b=polyder(p1,p2)p1*p2乘积的导数polyder(p1,p2)多项式值的计算polyval(p1,1.35)多点值的求法：a=3 4 5 6; polyval(p1,a)曲线拟合与插值一、 插值：1 一维线性插值插值定义为对数据点之间函数的估值方法，这些数据点是由某些集合给定。当人们不能很快地求出所需中间点的函数值时，插值是一个有价值的工具。最简单插值的例子是MATLAB的作图。按缺省，MATLAB用直线连接所用的数据点以作图。这个线性插值猜测中间值落在数据点之间的直线上。当然，当数据点个数的增加和它们之间距离的减小时，线性插值就更精确。MATLAB线性插值的直线可用函数interp1。yi=interp1(x, y, xi, methed )，这里x是独立变量(横坐标)，y是因变量(纵坐标)，xi是进行插值的一个数值数组,yi是对应的函数值method有以下几个参数：nearest:最近插值linear：线性插值spline：三次样条插值cubic：三次插值method的默认值为线性插值需要注意的是，对于上述的所有插值调用格式，人们不能要求有x变量范围以外的结果，而且要求向量x为单调，当x单调且等间距时为佳。例1x=linspace(0,2*pi,10);y=5,7,8,10,13,14,15,17,19,20;x1=1.2 2.1 3;y1=interp1(x,y,x1)2 Lagrang插值公式：Li(x)=(x-x0)(x-xi-1)(x-xi+1)(x-xn)(xi-x0)(xi-xi-1)(xi-xi+1)(x-xn)Ln(x)=yili(x)function s=lagrange(x,y,x0)%lagrange插值，x与y为已知的插值节点及其函数值，x0为需要求的插值点的x值。nx=length(x);m=length(x0);for i=1:m t=0.0; for j=1:nx u=1.0; for k=1:nx if k=j u=u*(x0(i)-x(k)/(x(j)-x(k); end end t=t+u*y(j); end s(m)=t;endreturn 见lagrange.m例2：对y=1/(1+x2),-5=x=5，用n=11个节点（等分）作上述三种插值，用m=21个插值点（等分）作图比较结果；x=linspace(-5,5,21);y=1./(1+x.2);x0=linspace(-5,5,11);y0=1./(1+x0.2);y1=lagrange(x0,y0,x);y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,spline);x,y,y1,y2,y3;plot(x,y,k:,x,y1,x,y2,x,y3)legend(true,lagrange,piece-liner,spline)gtext(lagr),gtext(piece_liner),gtext(spline),gtext(y=1/(x+x2) 见contrast.m3. 二维栅格点内插法：二维插值是基于与一维插值同样的基本思想。然而，正如名字所隐含的，二维插值是对两变量的函数z=f(x, y)进行插值。二维插值更为复杂，只是因为有更多的量要保持跟踪。interp2(x, y, z, xi, yi, method),这里x和y是两个独立变量，z是一个应变量矩阵。可选的参数method可以是:邻近点内插法（nearest）二维线性内插法（bilinear）二维样条内插法(spline)二维三次内插法(bicubic)x,y=meshgrid(-3:1:3);z=peaks(x,y);xi,yi=meshgrid(-3:0.25:3);zi1=interp2(x,y,z,xi,yi,nearest);zi2=interp2(x,y,z,xi,yi,bilinear);zi3=interp2(x,y,z,xi,yi,bicubic);subplot(2,2,1);surf(x,y,z);axis tight;title(original);% 先在三维坐标画出原始数据，看一下该数据的粗糙程度subplot(2,2,2);surf(xi,yi,zi1);axis tight;title(nearest);subplot(2,2,3);surf(xi,yi,zi2);axis tight;title(bilinear);subplot(2,2,4);surf(xi,yi,zi3);axis tight;title(bicubic);见shj.m绘制等高线（观察光滑度）x,y=meshgrid(-3:1:3);z=peaks(x,y);xi,yi=meshgrid(-3:0.25:3);zi1=interp2(x,y,x,xi,yi,nearest);zi2=interp2(x,y,z,xi,yi,bilinear);zi3=interp2(x,y,z,xi,yi,bicubic);subplot(2,2,1);contour(x,y,z,20);title(original);subplot(2,2,2);contour(xi,yi,zi1,20);title(nearest);subplot(2,2,3);contour(xi,yi,zi2,20);title(bilinear);subplot(2,2,4);contour(xi,yi,zi3,20);title(bicubic);4. 二位散布点内插法如果采样数据并非分布在栅格点上，我们利用zi=griddata(x,y,z,xi,yi)x=6*rand(100,1)-3;y=6*rand(100,1)-3;z=peaks(x,y);xi,yi=meshgrid(-3:0.2:3,-3:0.2:3); %用来产生三维绘图时的阵列（网格）zi=griddata(x,y,z,xi,yi);mesh(xi,yi,zi);hold onplot3(x,y,z,o);hold offaxis tight; %紧缩形式axis equal; %具有相等的刻度比例见dian.mview(2)表示正上方俯视图形每一圆球代表采样点（共100）而曲面则是使用griddata命令的内插结果，griddata命令只进行内插，而不能外插，因此所产生曲面数据只定义在采样点所形成的最小凸多边形内二、 数据拟合曲线拟合涉及回答两个基本问题：最佳拟合意味着什么？应该用什么样的曲线？可用许多不同的方法定义最佳拟合，并存在无穷数目的曲线。当最佳拟合被解释为在数据点的最小误差平方和，且所用的曲线限定为多项式时，那么曲线拟合是相当简捷的。数学上，称为多项式的最小二乘曲线拟合。对各数据点距离求平方，并把平方距离全加起来，就是误差平方和。这条虚线是使误差平方和尽可能小的曲线，即是最佳拟合。最小二乘这个术语仅仅是使误差平方和最小的省略说法。在MATLAB中，函数polyfit求解最小二乘曲线拟合问题。为了用polyfit，我们必须给函数赋予上面的数据和我们希望最佳拟合数据的多项式的阶次或度。如果我们选择n=1作为阶次，得到最简单的线性近似。通常称为线性回归。相反，如果我们选择n=2作为阶次，得到一个2阶多项式。现在，我们选择一个2阶多项式。(1) load census.matplot(cdate,pop,o) xlabel(年份)ylabel(人口（单位：百万）