电子科技大学《数学实验》2008-2009学年期末试题(含答案)

电子科技大学二零零八到二零零九学年第二学期期末考试《数学实验》课程考试题A卷(120分钟) 考试形式：闭卷考试日期：2009年7月8日
一、单项选择题(20分)
1、三阶幻方又称为九宫图，提取三阶幻方矩阵对角元并构造对角阵用( )
(A) diag(magic(3)); (B) diag(magic);
(C) diag(diag(magic(3))); (D) diag(diag(magic))。

2、MATLAB命令P=pascal(3)将创建三阶帕斯卡矩阵，max(P)的计算结果是( )
(A) 1 2 3 (B) 1 2 1 (C) 3 6 10 (D) 1 3 6
3、命令J=*1;1;1+**1,2,3+;A=j+j’-1将创建矩阵( )
(A)
123
234
345
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎣⎦
; (B)
234
345
456
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎣⎦
(C)
123
123
123
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎣⎦
(D)
111
222
333
⎡⎤
⎢⎥
⎢⎥
⎢⎥
⎣⎦
4、data=rand(1000,2);x=data(:,1);y=data(:,2);II=find(y<sqrt(x)&y>x.^2);的功能是( )
(A) 统计2000个随机点中落入特殊区域的点的索引值；
(B) 统计1000个随机点落入特殊区域的点的索引值；
(C) 模拟2000个随机点落入特殊区域的过程；
(D) 模拟1000个随机点落入特殊区域的过程。

5、MATLAB计算二项分布随机变量分布律的方法是( )
(A) binocdf(x,n,p); (B) normpdf(x,mu,s); (C)binopdf(x,n,p); (D) binornd(x,n,p)。

6、MATLAB命令syms e2;f=sqrt(1-e2*cos(t)^2);S=int(f,t,0,pi/2)功能是()
(A) 计算f(x)在[0,pi/2]上的积分；(B) 计算f(t)不定积分符号结果；
(C) 计算f(x)积分的数值结果；(D) 计算f(t)定积分的符号结果。

7、y=dsolve(‘Dy=1/(1+x^2)-2*y^2’,’y(0)=0’,’x’);ezplot(y)的功能是( )
(A) 求微分方程特解并绘图；(B) 解代数方程(C) 求定积分；(D)求微分方程通解。

8、X=10000 ;0.5*asin(9.8*X/(515^2))的功能是计算关于抛射体问题的()
(A) 十公里发射角；(B) 十公里飞行时间；(C)最大飞行时间；(D)最大射程。

9、theta=linspace(0,2*pi,100) ;r=cos(4*theta) ;polar(theta,r,’k’)功能是()
(A) 绘四叶玫瑰线；(B)绘三叶玫瑰线；(C)绘心脏线；(D) 绘八叶玫瑰线。

10、北京和纽约的经度分别是：东经118和西经76，根据经度差计算时差用()
(A) fai1=118;fai2=-76;Dfai=(fai1+fai2)/24; (B) fai1=118;fai2=-76;Dfai=(fai1+fai2)/15;
(C) fai1=118;fai2=-76;Dfai=(fai1-fai2)/24; (D) fai1=118;fai2=-76;Dfai=(fai1-fai2)/15。

二、程序阅读题(40分)
1、直方图功能是将数据分为n个类，统计各个类的数据量并绘图。

借用现有的直方图命令hist，编写新直方图程序如下。

function m=myhist(data,n)
if nargin==1,n=7;end
Xmin=min(data);Xmax=max(data);h=(Xmax-Xmin)/n;
m=hist(data,n)/length(data)/h;
t=linspace(Xmin,Xmax,n+1);
II=1:4:4*n-3;JJ=1:n;
x(II)=t(JJ);y(II)=zeros(1,n);
x(II+1)=t(JJ);y(II+1)=m;
x(II+2)=t(JJ+1);y(II+2)=m;
x(II+3)=t(JJ+1);y(II+3)=zeros(1,n);
plot(x,y,'k')
(1) 变量data存放了1000个数据，在命令窗口调用myhist(data)的结果是( )
(A) 只绘数据的直方图而不显示被分类后各类的数据量；
(B) 只显示被分类后各类的数据量而不绘数据的直方图；
(C) 既绘数据直方图也显示被分类后各类的数据量；
(D) 根据默认值在数据范围内插入七等分点绘直方图。

(2) 关于新直方图绘图程序下面说法不正确的是( )
(A)h是n等分直方图中小区间长度；(B) 修改程序最后一行可绘红色直方图；
(C) 直方图中所有小矩形面积之和为1；(D) 直方图中所有小矩形的高度和为1。

2、3n+1问题反映一个数学猜想：对任一自然数n，按如下法则进行运算：若n为偶数，则将n除2，若n为奇数，则将n乘3加1。

重复这种操作，结果终会为1。

实验程序如下。

function [k,N]=threeN(n)
if nargin==0,n=5;end
k=1;N=n;
while n~=1
r=rem(n,2);
if r==0
n=n/2;
else
n=3*n+1;
end
N=[N,n];k=k+1;
end
(1)在MATLAB命令窗口中直接调用threeN运行结果为( ) %5 16 8 4 2 1
(A)只显示k的最后数值为6；(B) 只显示k的最后数值5；
(C) 同时显示k和N的数据；(D) 仅显示N的所有数据。

(2)实验程序运行过程中( )
(A) 输入变量n不发生改变；(B)N是记录数据变化的一维数组；
(C) N记录每次数据变化的单个数据；(D)n是记录数据变化的一维数组。

3、将半径为r的球体(密度1
ρ)置入水中，球体将浮出水面一定高度。

程序如下：
<
function [h,Rou]=highNu(r)
if nargin==0,r=10;end
Rou=0.3:0.1:1;
N=length(Rou);
for k=1:N
rouk=Rou(k);
P=[1,-3*r,0,4*r^3*rouk];
x=roots(P);
II=find(x<2*r&x>0); h(k)=2*r-x(II); end
(1)在MATLAB 命令窗口省略输入调用函数highNu 将显示( )
(A) 球体浮出水面的高度数据； (B) 球体的8个不同的密度数据； (C) 球体沉入水下的深度数据； (D) 深度数据和密度数据。

(2) 程序中变量x 存入如下方程的根( )
(A)043323=+-ρr rx x ； (B)ρ3343r r x =-； (C)033=-r x ； (D)0313=-rx 4、一阶常微分方程确定一个平面向量场，初值条件确定了向量场中一条曲线。

程序如下：
[x,y]=meshgrid(0:.25:6,0:.05:2); k=y.*(1-y); d=sqrt(1+k.^2); px=1./d;py=k./d; quiver(x,y,px,py),hold on u=dsolve('Du=u*(1-u)','u(0)=.2'); v=dsolve('Dv=v*(1-v)','v(0)=1.8'); ezplot(u,[0,6]) ezplot(v,[0,6])
(1) 程序中所绘向量场对应的一阶常微分方程是( ) (A)2
1y y '+=
'； (B))1(y y y -='； (C)2
1/1y y '+='； (D))1(/1y y y -='。

(2) 关于实验程序下面说法错误的是( )
(A) 程序中第一个初值条件所对应的解曲线在图1中上方； (B) 程序中第二个初值条件所对应的解曲线在图1中上方； (C) 程序绘图原理是根据每一点处曲线切线的单位向量绘图；
(D) 当初值数据大于1时解曲线单调减少，当初值数据小于1时解曲线单调增加。

5、维维安尼体由柱面切割球体所得。

下面程序的功能是演示柱面切割球体的过程。

function viviani(dt) if nargin==0,dt=10;end N=fix(360/dt); [X,Y,Z]=sphere(N); mesh(X,Y,Z),hold on [x,y,z]=cylinder([1,1],N);
y=.5*y;x=.5*(1-x);z(1,:)=-ones(1,N+1); for p=10:N+1 II=1:p;u=x(:,II); v=y(:,II);w=z(:,II); mesh(u,v,w),pause(.5) end
(1) 根据程序中语句，所绘图形中( )
(A) 圆柱的半径为1； (B) 圆柱的高度为1； (C) 圆柱以Z 轴对称； (D) 圆柱的高度为2。

(2) 关于实验程序以下错误的说法是( )
t
-9/(-9+4 exp(-t))
图1 向量场图
(A) 程序中输入变量dt大则球面网格线稀；(B) 程序正常运行时球面图形保持不变；
(C) 程序绘图时每半秒种图形变动一次；(D) 每循环一次只加绘柱面一条母线。

三、程序填空(40分)
1、中国农历60年一大轮回，按天干“甲乙丙丁戊已庚辛壬癸”和地支“子丑寅卯辰巳午未申酉戍亥”循环排列而成。

已知2009年是农历已丑年，通过简单计算可以找出年份与天干/地支对应的规律。

下面数学实验程序对输入年份，计算并输出字符串农历纪年。

填空完善程序。

function calendar=year(year)
ifnargin==0, year=2009;end
S1=’ 甲乙丙丁戊已庚辛壬癸’;
S2=’子丑寅卯辰巳午未申酉戍亥’;
k1=①; %定位天干序数%2010庚寅
s1=S1(k1);
k2=②; %定位地支序数
s2=S2(k2);
calendar=strcat(int2str(year),’年是’,s1,s2,’年’)
2、红、绿两队从相距100公里的地点同时出发相向行军。

红队速度为10(公里/小时)，绿队速度为8(公里/小时)。

开始时，通讯员骑摩托从红队出发为行进中的两队传递消息。

摩托车的速度为60(公里/小时)往返于两队之间。

每遇一队，立即回驶向另一队。

当两队距离小于0.2公里时，摩托车停止，下面数学实验程序模拟计算摩托车跑了多少趟。

请填空完善程序。

function k=moto(A,B)
ifnargin==0,A=0;B=100;end
va=10;vb=8;vc=60;
f=1;k=0;
while (B-A)>0.2
if f==1
tk=(B-A)/(vb+vc);
else
tk=①; %计算A与C相遇时间
end
A=②; %计算A点位置
B=③; %计算B点位置
f=-f;
k=k+1;
end
3、为了进入地月转移轨道，嫦娥一号卫星进行了四次变轨调速度。

第一次变轨从16小时初始轨道进入16小时轨道，第二次卫星进入24小时轨道，第三次卫星进入48小时轨道，第四次卫星进入116小时地月转移轨道。

上面小时数并不是准确轨道周期，变轨目的是将速度从10.3(km/s)逐渐提高到约10.9(km/s)。

下面数学实验程序是在区间[10.3,10.9]上插入线性等分点，即每个轨道的最大速度以等差数列出现，然后近似计算出每个轨道的周期参数。

填空完善程序。

function satel1()
R=6378;
h=[200,600,600,600,600];
H=[51000,51000,71000,128000,370000];
a=(h+H+2*R)/2;
c=(H-h)/2;
b=①; %计算短半轴数据
E2=(c./a).^2;
L=2*pi*a.*(1-E2/4-3*E2.^2/64)
format bank
Vmax=linspace(10.3,10.9,5)
S=②; %根据最大速度计算每秒钟扫过的面积
Times=a.*b.*pi./S;
myTimes=Times/3600
4、冰淇淋锥的下部为圆锥面，上部为半球面。

计算体积的蒙特卡罗方法是在包含冰淇淋的六面体内产生N个均匀分布的随机点，并统计落入锥体内的随机点的数目m。

根据比值m/N 和六面体体积数据计算出锥体体积数据，这种随机统计方法会产生误差，根据大数定律，误差变量服从正态分布。

下面数学实验程序使用上面二题中第1小题绘出误差直方图与正太分布密度函数比较，填空完善程序。

if nargin==0,L=1000;end
for k=1:L
P=rand(2000,3);
x=2*P(:,1)-1;
y=2*P(:,2)-1;
z=①; %计算随机点Z坐标数据
R2=x.^2+y.^2;
R=) ②; %计算随机点到坐标原点距离II=find(z>=R&z<=1+sqrt(1-R2));
m=length(II);
q(k)=8*m/2000;
end
X=q-pi;
mu=mean(X);
sagma=sqrt(sum((X-mu).^2)/(L-1));
myhist(X,7);hold on
x=linspace(-3*sagma,3*sagma,50);
y=③; %计算正态分布密度函数值
plot(x,y,'r')
附参考答案：
一、单项选择题(每小题2分共20分)
CDAB CDAA DD
二、程序阅读题(每小题4分共40分)
1、CD
2、AB
3、AA
4、BB
5、DD
三、程序填空(每小题4分共40分)
1、mod(year-4,10)+1;
mod(year-4,12)+1; 2、(B-A)/(vc+va); A+va*tk; B-vb*tk;
3、sqrt(a.*a-c.*c); (R+h).*Vmax/2;
4、2*P(:,3);
sqrt(R2);
normpdf(x,mu,sagma) 或 exp(-(x-mu).^2/(2*sagma^2))/(sqrt(2*pi)*sagma) 22
2)(21σ
μσ
π--x e
电子科技大学二零零七到二零零八学年第二学期期末考试
《数学实验》课程考试题 A 卷 (120分钟) 考试形式：闭卷考试日期：2008年6月27日
所有答案一律写在答题纸上，写在试卷上无效。

一、单项选择题(20分)
1．MA TLAB 命令A=rand(5,5)；创建55()ij A a ´=，求5
1
m a x
||ij j
i a =å
用()
(A) max(sum(abs(A))); (B) max(sum(abs(A ’))); (C) max(sum(A))); (D) sum(max(A)); 2．MA TLAB 命令x=[1,2,4,5,9];mean(x)，的计算结果是( ) (A) 4 (B) 4.2 (B) 4.5 (D) 21 3．MA TLAB 命令x=rand(10,1)生成10个随机数，将它们从大到小排序，使用( ) (A) y=sort(x);z=y(10:1); (B) [y,II]=sort(x);z=y(II); (C) y=sort(x);z=y(10:-1;1); (D) [y ,II]=sort(x);z=x(II); 4．MA TLAB 命令roots([1,0,0,-1])的功能是() (A) 产生向量[1,0,0,1]； (B) 求方程310x +=的根； (C) 求多项式31x -的值 (D) 求方程310x -=的根。

5．MA TLAB 命令A=magic(3)创建3阶幻方矩阵，求A 的特征值绝对值最小用( ) (A) min(abs(eig(A))); (B) min(eig(abs(A))); (C)min(eig(A)); (D) min(abs(A)); 6．命令factor()用于分解因式，syms x; f=4*x^3+9*x^2-30*x; factor(diff(f))的结果是() (A) (x-1)*(2*x-5) (B) 6*(x-1)*(2*x+5) (C) 6*(x+1)*(2*x+5) (D) (x+1)*(2*x-5) 7．MA TLAB 命令syms x; f=sin(x); V=pi*int(f*f,x,0,pi)功能是() (A) 绘出函数f 在[0,2p ]图形； (B) 计算函数f 在[0,2p ]的积分； (C) 计算旋转曲面所围的体积； (D) 计算旋转曲面的表面积。

8．十二属相为“鼠牛虎兔龙蛇马羊猴鸡狗猪”，命令k=mod(2008,12)+1的结果是( )
(A) k指向第二动物牛；(B) k指向第三动物虎；(C) k指向第四动物兔；(D) k指向第五动物龙。

9．MA TLAB命令[x,y]=meshgrid(1:3);H=1./(x+y-1)产生的矩阵H是( )
(A)
111
222
333
轾
犏
犏
犏
犏
犏臌
(B)
123
123
123
轾
犏
犏
犏
犏
犏臌
(C)
123
233
345
轾
犏
犏
犏
犏
犏臌
(D)
11/21/3
1/21/31/4
1/31/41/5
轾
犏
犏
犏
犏
犏
臌
10．下面有关MA TLAB变量名和函数名的说法中，错误的说法是( )
(A) 变量名的第一个字符必须是一个英文字母
(B) 变量名可由英文字母、数字和下划线混合组成
(C) 变量名不得包含空格和标点，但可以有下连字符
(D) 变量名和函数名对于英文的大小使用没有区别
二、程序阅读题(40分)
1．传说古希腊曾流行瘟疫，人们为消除灾难求助于神。

神说：把神庙中黄金祭台增容一倍，可消除瘟疫。

当立方体祭台尺寸放大一倍后，瘟疫仍然流行。

人们才知道体积并不是扩大了两倍。

这个古希腊难题被称为倍立方体问题，在人类还没有认识到无理数时，企业界企图用
文本框所示
a=2^(1/3);
D=1;
for k=1:8
D=D*10;
b=fix(a*D)/D;
V(k)=b^3;
end
error=V’-2
(1) 程序中循环控制变量k从1变量8，而变量D=10k的作用是( )
(A) 将a的小数点向右移D位取整；
(B) 将a的小数点向右移D位取整后再向左移D位；
(C) 将a的小数点向右移k位取整后再向左移k位；
(D) 将a的小数点向左移k位取整后再向右移k位；
(2) 程序中变量b存放的数据是( )
(A) 将a的小数点后第k位减1所得；(B) 将a的小数点k位后按四舍五入所得；
(C) 将a的小数点后第k位增1所得；(D) 将a的小数点k位后截断舍去所得。

2．Viviani体是圆柱体222
(/2)/4
x R y R
-+被球面2222
x y z R
++=所割立体。

下面的数学实验程序功能是取R=2求体积上半部分，先利用符号处理重积分并转换为数值数据，再用蒙特卡罗方法计算体做对比。

数学实验程序如下：
syms x y;
f=sqrt(4-x^2-y^2);
y1=sqrt(2*x-x^2);y2=sqrt(2*x-x^2);
S1=int(f,y,y1,y2);S2=int(S1,x,0,2)
V=double(S2)
图1 V ivinai问题
P=rand(10000,3);
X=2*P(:,1);Y=2*P(:,2);Z=2*P(:,3); II=find((X-1).^2+Y .^2<=1&Z<=sqrt(4-X.^2-Y .^2)); V1=8*length(II)/10000
(1) 符号计算所用的积分公式是( )
(A) 2
V d x
y -=
蝌 (B) 210
1V dy
-
=
蝌
(C) V y -
=
ò
(D) V x -
=
ò
(2) 蒙特卡罗方法选用的随机点变化范围的立方体区域是( ) (A) {(,,)|(0,2),(0,2),(0,2)}x y z x y z W =挝 ； (B) {(,,)|(0,2),(1,1),(0,2)}x y z x y
z W =挝-
(C) {(,,)|02),01,02}x y z x y z W =<<<<<< (D) {(,,)|02),02,02}x y z x y z W =<<<<<<
3．某厂生产两种产品，产一吨甲产品用A 资源3吨、B 资源4m 3
；产一吨乙产品用A 资源2吨，B 资源6m 3，C 资源7个单位。

一吨甲产品和乙产品分别价值7万元和5万元，三种资源限制分别为90吨、200m 3
和210个单位。

生产两种产品使总价值最高的生产方案可用数学实验程序计算。

C=[-7,-5];A=[3 2;4 6;0 7];b=[90;200;210]; Aeq=[];Beq=[];
e0=[0,0];e1=[inf,inf];
[x,fval]=linprog(C,A,b,Aeq,beq,e0,e1); (1) 程序中变量C 表示( )
(A) 目标函数系数； (B) 等式约束系数； (C) 不等式约束系数； (D) 等式约束常向量 (2) 程序中变量A 表示( )
(A) 等式约束矩阵； (B) 不等式约束矩阵； (C) 决策变量的值； (D) 目标函数的最大值 4．用十二星座反映人的心理和行为。

十二星座是：白羊座、金牛座、双子座、巨蟹座、狮子座、处女座、天秤座、天蝎座、射手座、魔蝎座、水瓶座、双鱼座。

游戏规则如下：确定一个正整数k(0<k<13)对应星座之一，将四颗骰子同时掷一次，由点数之和确定游戏者是否是第k 个星座。

模拟程序如下： function Fn=playingstar(k)
if nargin==0,k=2;end
S='白羊座金牛座双子座巨蟹座狮子座处女座天秤座天蝎座射手座魔蝎座水瓶座双鱼座'; if k<1|k>12,error('please input again 1 to 12');end k1=3*(k-1)+1;k2=3*k; Sk=S(k1:k2)
Show=strcat('你选择了----',Sk) N=2000;R=1+fix(6*rand(4,N)); x=sum(R);y=mod(x,12)+1;
II=find(y==k); %第十行语句
n=length(II);Fn=n/N
(1) 当用户调用函数程序时，没有输入数据，则程序运行后将显示2000次随机实验( ) (A) 游戏者可能是白羊座的频率； (B) 游戏者可能是金牛座的频率； (C) 游戏者可能是双子座的频率； (D) 游戏者可能是其它星座的频率。

(2) 第十行语句的功能是( )
(A) 统计2000次随机实验中，游戏者可能是第k 个星座的频率； (B) 统计2000次随机实验中，游戏者可能是第k 个星座的索引值； (C) 统计2000次随机实验中，游戏者可能是第k 个星座的次数； (D) 统计2000次随机实验中，游戏者可能是第k 个星座的频数。

5．一个平面多边形由它的n 个顶点确定，将顶点按逆时针方向排列为：111(,)P x y ，...，
(,)n n n P x y 。

将第(n+1)个顶点设为111(,)P x y 。

则多边形面积可由二阶行列式求和计算，数
学实验程序如下：
data=[-1,-1;1,-1;1,1;0,0;-1,1]; n=size(data,1);Sk=0;pk1=data(1,:); for k=2:n
pk=data(k,:);Dk=det([pk1;pk]); Sk=Sk+Dk;pk1=pk; end
pk=data(1,:);Dk=det([pk1;pk]); Sk=Sk+Dk; Sn=0.5*Sk
(1) 程序中所用的二阶行列式是( ) (A) 11k k k k k x y D x y ++=
；(B)1
1
k k k k k x y D x y --=
；
(C)1
1
k k k k
k
x y D x y --=；(D)11k k k k
k
x y D x y ++=
(2) 程序中所用的多边形求和公式是( ) (A) 1
1
2
n n k k S D ==
å；(B)1
n
n k k S D ==
å
；(C)12
1
2
n n k k S D +==
å；(D)1
2
n n k k S D +==
å
三、程序填空(10分)
1．二阶正交矩阵作用于某一向量时，其效果是将该向量旋转，旋转解为a (逆时针旋转为正)。

把一个以原点为中心的正三角形旋转/50p ，并缩小90%，迭代33次创建图3。

完成程序填空：
bata=[1/2;7/6;11/6;15/6]*pi; x=cos(bata);y=sin(bata); line(x,y) xy=[x,y];
图2 多边形面积计算
alfa=pi/50;
A=[cos(alfa),-sin(alfa);sin(alfa),cos(alfa)]; for k=1:33 xy=①; x=xy(:,1); y=②; line(x,y) end
2．长征三号甲运载火箭提供给探月卫星的初始速度不足以将卫星送往月球轨道。

为提高到奔月速度，中国航天工程师使用了卫星变轨技术。

数学实验程序根据变轨中轨道周期和近地点距离数据，利用开普列第二定律模拟计算计算卫星飞行的最大速度。

填空完善下面实验程序。

R=6378;Time=[16,15.63,23.3,50.5,225]*3600;
h=[200,600,600,600,600];H=[51000,51000,71000,128000,370000]; a=(h+H+2*R)/2; c=①
b=sqrt(a.*a-c.*c); S=②
Vmax=2*S./(R+h)
3．抛射体运动的伽利略模型是针对无阻力情况。

考虑阻力与速度成正比的数学模型，可导出参数方程
002
sin cos (1ex p ())/ (
)(1ex p ())v g g x v kt k y kt t k
k
k
a
a =--=+
---，
完成程序填空：
k=0.1; alpha=pi/4; v0=198;g=9.8; t=0;dt=0.1;x=0;y=0; while y>=0 t=t+dt; xk=①
yk=((v0*sin(alpha)+g/k)*(1-exp(-k*t))-g*t)/k; x=[x,xk]; y=② end plot(x,y,'r')
4．五月十二日以来汶川地区发生五级以上地震已经越过了25次，将地震数据整理，表示成n 行三列的矩阵data ，第二列为经度k x ，第一列为为纬度k y ，第三列为震级k d 。

不同震发点和震级的地震都对周边地区产生影响，距离近则影响强烈，距离远则影响减弱。

用地震影响曲面描述，其数学原理如下：以每次震发点的经纬度为中心构造函数 22
()()(,)exp (), (1,2,....,)k k k k k
x x y y z x y d k n d -+-=-=
图4 有阻力的抛射曲线
将这一函数离散化为矩阵并逐次累加，最后除以累加后的矩阵的最大值，可绘出地震影响曲面如图5。

数学实验程序如下，请填空完善。

load data.txt
d=data(:,3);n=length(d);
x=data(:,2);
y=data(:,1);
[X,Y]=meshgrid(100:0.2:110,30:0.2:35); Z=zeros(size(X));
for k=1:n
xk=x(k);yk=y(k);
dk=d(k);
Z=Z+①
end
Maxz=②
Z=Z./Maxz;
mesh(X,Y,Z)
colormap([0 0 0])
图5 地震影响曲面。