数学实验与Matlab
【 clf,x=linspace(-2,2,40); y=x; z=x;
[X,Y,Z]=meshgrid(x,y,z); w=X.^2+Y.^2+Z.^2;
slice(X,Y,Z,w,[-1,0,1],[-1,0,1],[-1,0,1]),colorbar 】
11
3. 空间曲线及其运动方向的表现：plot3和quiver指令
【 clf, t=0:0.1:1.5;
Vx=2*t;Vy=2*t.^2;Vz=6*t.^3-t.^2;
x=t.^2;y=(2/3)*t.^3;z=(6/4)*t.^4-(1/3)*t.^3; %由速度得到曲线
plot3(x,y,z,'r.-'),hold on %画飞行
轨迹
%算数值梯度，也就是重新计算数值速度矢量，这只是为了编程的方便，不是必
须的
Vx=gradient(x);Vy=gradient(y);Vz=gradient(z);
quiver3(x,y,z,Vx,Vy,Vz),grid on %画速度矢量图 xlabel('x'),ylabel('y'),zlabel('z') 】
5z102.2应用、思考和练习
2.2.1 线性p周期函数
zxy2_3_f.m
【 function f=zxy2_3_f(x)
042y001x32 图2.12 飞机的飞行轨迹与方向
f=sin(x+cos(x)); 】 zxy2_3.m 【 clear，clf
a=-8;b=12;n=300;xx=linspace(a,b,n); h=zxy2_3_f(xx); S(1)=0; for i=2:n
S(i)=S(i-1)+quad('zxy2_3_f',xx(i-1),xx(i));
11
数学实验与Matlab
end
subplot(1,2,1),plot(xx,S,'k-'),axis([a,b,-1.5,9])
12
subplot(1,2,2),plot(xx,[h;zeros(1,length(xx))],'k-'),axis([a,b,-1.5,1.5]) 】
2.2.2 平面截割法和曲面交线的绘制
◆用平行截面法讨论由曲面 z?x2?2y2 构成的马鞍面形状。
zxy2_6.m （ 平行截割法示例 ， 本程序的绘制两曲面交线方法可以套用） 【 clf, a=-20;eps0=1;
[x,y]=meshgrid(-10:0.2:10); %生成平面网格 v=[-10 10 -10 10 -100 100]; %设定空间坐标系的范围
colormap(gray) %将当前的颜色设置为灰色 z1=(x.^2-2*y.^2)+eps; %计算马鞍面函数z1=z1(x,y) z2=a*ones(size(x)); %计算平面 z2=z2(x,y) r0=abs(z1-z2)<=eps0;
%计算一个和z1同维的函数r0,当abs(z1-z2)<=eps时r0 =1；当abs(z1-z2)>eps0时，r0 =0。
%可用mesh(x,y,r0)语句观察它的图形，体会它的作用，该方法可以套用。
zz=r0.*z2;xx=r0.*x;yy=r0.*y; %计算截割的双曲线及其对应的坐标 subplot(2,2,2), %在第2图形窗口绘制双曲线
h1=plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'+'); set(h1,'markersize',2),hold on,axis(v),grid on subplot(2,2,1), %在第一图形窗口绘制马鞍面和平面 mesh(x,y,z1);grid,hold on;mesh(x,y,z2); h2=plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'.'); %画出二者的交线 set(h2,'markersize',6),hold on,axis(v),
for i=1:5 %以下程序和上面是类似的，通过循环绘制一系列的平面
去截割马鞍面
a=70-i*30; %在这里改变截割平面
z2=a*ones(size(x)); r0=abs(z1-z2)<=1; zz=r0.*z2;yy=r0.*y;xx=r0.*x; subplot(2,2,3),
mesh(x,y,z1);grid,hold on;mesh(x,y,z2);hidden off
12
数学实验与Matlab
h2=plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'.'); axis(v),grid subplot(2,2,4),
h4=plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'o'); set(h4,'markersize',2),hold on,axis(v),grid on end 】
13
2.2.3 微分方程的斜率场
◆ 绘制微分方程 场和解曲线的关系。
zxy2_5.m ( 绘制一阶微分方程的斜率场和解曲线)
【 clf,clear %清除当前所有图形窗口的图像，清除当前工作空间的
内存变量。
dy并将解曲线画在图中，观察斜率?xy, y(0)?0.4的斜率场，
dta=0;b=4;c=0;d=4;n=15;
[X,Y]=meshgrid(linspace(a,b,n),linspace(c,d,n)); %生成区域中的网格。 z=X.*Y;
%计算斜率函数。
Fx=cos(atan(X.*Y));Fy=sqrt(1-Fx.^2); %计算切线斜率矢量。 quiver(X,Y,Fx,Fy,0.5),hold on,axis([a,b,c,d])
%在每一网格点画出相应的斜率矢量，0.5是控制矢量大小的控制参数，可以调整。
[x,y]=ode45('zxy2_5f',[0,4],0.4); %求解微分方程。
%zxy2_5f.m是方程相应函数f(x,y)的程序，单独编制；[x0,xs]=[0,4]为求解区间；
%y0=0.4为初始值；输出变量x,y分别为解轨线自变量和因变量数组。
plot(x,y,'r.-') %画解轨线 】 zxy2_5f.m （微分方程的函数子程序） 【 function dy=zxy2_5f(x,y)
dy=x.*y; 】
13
数学实验与Matlab
14
2.2.4颜色控制和渲染及特殊绘图指令
1.地球表面的气温分布（sphere指令）
◆
【 [a,b,c]=sphere(40);t=max(max(abs(c)))-abs(c);surf(a,b,c,t);
axis('equal'),colormap('hot'), shading flat,colorbar 】
2.旋转曲面的生成：柱面指令cylinder和光照控制指令surfl
◆
【 x=0:0.1:10;z=x;y=1./(x.^3-2.*x+4);
[u,v,w]=cylinder(y);surfl(u,v,w,[45,45]); shading interp 】
3.若干特殊图形
◆ 运行下面程序，了解各指令的用法和效果。 【 x=[1:10]; y=[5 6 3 4 8 1 10 3 5 6];
subplot(2,3,1),bar(x,y),axis([1 10 1 11]) subplot(2,3,2),hist(y,x),axis([1 10 1 4]) subplot(2,3,3),stem(x,y,'k'),axis([1 10 1 11]) subplot(2,3,4),stairs(x,y,'k'), axis([1 10 1 11])
subplot(2,3,5), x = [1 3 0.5 5];explode = [0 0 0 1];pie(x,explode) subplot(2,3,6),z=0:0.1:100; x=sin(z);y=cos(z).*10; comet3(x,y,z) 】
实验３．函数式－直接确定型模型
3.1知识要点和背景：函数 — 直接确定性模型
14
数学实验与Matlab
15
3.2实验与观察：插值与拟合
3.2.1 插值方法与多项式拟合的概念 3.2.2 用Matlab作插值和拟合
3.2.3 用鼠标选节点 观察插值、拟合的效果
3.2.4 观察程序说明
zxy3_1.m
【 clf,a=-1;b=1;n=100;
%用内联函数inline命令定义函数，
%在后面可直接用于函数g的计算,要改变函数做实验，可按此格式重新定义g
g=inline('x^2-x^4'); xx=linspace(a,b,n); for i=1:n
gx(i)=gxx(i)); % 前面已经用inline命令定义了g，可以这样用g计算函
数值
end
ymin=min(gx)*0.8;ymax=max(gx)*1.2;%分四个界面画图g的图形，以便于结
果比较
subplot(2,2,1),
plot(xx,gx,'--'),grid,hold on,axis([a b ymin ymax]),title('近邻插值') subplot(2,2,2),
plot(xx,gx,'--'),grid,hold on,axis([a b ymin ymax]),title('线性插值') subplot(2,2,3),
plot(xx,gx,'--'),grid,hold on,axis([a b ymin ymax]),title('样条插值') subplot(2,2,4),plot(xx,gx,'--'),
grid,hold on,axis([a b ymin ymax]),title('多项式拟合')
%用鼠标在屏幕上选点[x,y,button] = ginput(n)，可套用下面程序的格
15