中南大学
MATLAB程序设计实践
班级：材料1207 学号：********** 姓名：*** 时间：2014-6-26
基础题
（1）用Slice函数给出其整体分布特征；
A.流程图:
开始 从Data.txt中读取数据 通过循环形成三维矩阵 通过meshgrid函数生成绘制3-D图形所需的网格数据 用通过Slice函数给出其整体分布特征 结束 B.代码:
fid=fopen('Data.txt'); head=fscanf(fid,'%s',40);
for i=1:19 ch=fscanf(fid,'%f',1); d=fscanf(fid,'%f',[19,19]); data(:,:,i)=d'; end
[phi1,phi,phi2]=meshgrid(0:5:90,0:5:90,0:5:90); slice(phi1,phi,phi2,data,[45,90],[45,90],[0,45])
C.结果：
100806040200100500020406080100（2）用pcolor或contour函数分别给出(φ2＝0, 5, 10, 15, 20, 25, 30, 35 … 90)切面上f分布情况(需要用到subplot函数)；
A.流程图:
开始 从Data.txt中读取数据 通过meshgrid函数绘制二维网格 通过subplot函数生成19个切面 用通过contour函数给出各个切面上f分布情况 结束
B.代码:
fid=fopen('Data.txt'); head=fscanf(fid,'%s',40);
[phi1,phi]=meshgrid(0:5:90,0:5:90); for i=1:19
ch=fscanf(fid,'%f',1);
da=fscanf(fid,'%f',[19,19]); subplot(5,4,i) contour(phi1,phi,da') end
C.结果：
505050500050005000500050505050500050005000500050505050500050005000500050505050500050005000500050505050005000500050
(3) 用plot函数给出沿α取向线(φ1=0~90，φ＝45，φ2＝分布情况。
A.流程图:
开始 从Data.txt中读取数据 生成三维矩阵 plot函数给出沿α取向线(φ1=0~90，φ＝45，φ2＝0)的f分布情况。 结束 0)的f
B.代码:
fid=fopen('Data.txt'); head=fscanf(fid,'%s',40);
for i=1:19 d=fscanf(fid,'%f',[19,19]); data(:,:,i)=d'; end
b=data(10,:,1) a=0:5:90; plot(a,b,'-o')
C.结果：
1086420-20102030405060708090
科学计算部分
一、编程实现以下科学计算算法，并举一例应用之。
“矩形法、梯形法数值积分”
1.矩形法数值积分
A.算法说明：
矩形法是一种计算定积分近似值的方法，其思想是求若干个矩形的面积之和，这些矩形的高由函数值来决定。将积分区间[a, b] 划分为n个长度相等的子区间，每个子区间的长度为(a-b)/n 。这些矩形左上角、右上角或顶边中点在被积函数上。这样，这些矩形的面积之和就约等于定积分的近似值。
由函数上的点为矩形的左上角、右上角或顶边中点来决定，又分别被称为下（左）矩形公式、上（右）矩形公式和中矩形公式。当 n 逐渐扩大时，此近似值更加准确。上述的点无论取哪个值，最终和式的值都将趋近于定积分的值。 B.流程图：
开始 n=1;h=b-a 此时划分为1个子区间 I1=0;I2= (subs (sym (f), findsym (sym (f)), a) + subs (sym (f), findsym (sym (f) ) , b)) / h; 是 否 abs 否 n=n+1; h=(b-a) / n; I1=I2; for i=0:n-1 是 x=a+h*i; x1=x+h; I2=I2+h*subs(sym(f),findsym(sym(f)),x1); I=I2;Step=n; 结束 C.代码
function[I,step]=CombineTraprl(f,a,b,eps) % 复合矩形公式求函数f在区间[a,b]上的定积分 %函数名：f %积分下限：a %积分上限：b %积分精度：eps %积分值：I
%积分划分的子区间个数：step if(nargin==3) %如果参数个数为3 eps=1.0e-4; %默认精度为0.0001 end n=1; h=b-a; I1=0;
I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>eps n=n+1 h=(b-a)/n; I1=I2; I2=0;
for i=0:n-1 %第年n次的复合矩形公式积分
x=a+h*i; %i=0 和n-1时，分别代表积分区间的左右端点 x1=x+h
I2=I2+h*subs(sym(f),findsym(sym(f)),x1); end end I=I2; step=n;
D.应用举列：
复合矩形法求数值积分应用举例，利用复合矩形法计算定积分x2dx
12a.程序代码：
[q,s]=CombineTraprl('x^2',1,2) %精度为默认的10-4 b.结果：
2. 梯形法数值积分
A．算法说明：
梯形法数值积分采用的梯形公式是最简单的数值积分公式，函数f(x)在区间
[a,b]上计算梯形法数值积分表达式为：
bba f(x)dx[f(a)f(b)]
a2由于用梯形公式来求积分十分粗糙，误差也比较大，后来改进后提出了复合梯形公式：
hban,其中，n为积分区间划分的个数；h为积分步长。
在MATLAB中编程实现的复合梯形公式的函数为：Combine Traprl. 功能：复合梯形公式求函数的数值积分。
调用格式：[I,step]=CombineTraprl(f,a,b,eps). 其中，f为函数名； a为积分下限； b为积分上限； eps为积分精度； I为积分值；
Step为积分划分的区间个数 B.代码
function[I,step]=CombineTraprl(f,a,b,eps) % 复合梯形公式求函数f在区间[a,b]上的定积分 %函数名：f %积分下限：a %积分上限：b %积分精度：eps %积分值：I
%积分划分的子区间个数：step if(nargin==3)
eps=1.0e-4; %默认精度为0.0001 end n=1; h=(b-a)/2; I1=0;
I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; while abs(I2-I1)>eps n=n+1 h=(b-a)/n; I1=I2; I2=0;
for i=0:n-1 %第年n次的复合梯形公式积分
x=a+h*i; %i=0 和n-1时，分别代表积分区间的左右端点
x1=x+h
I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1)); end end I=I2; Step=n;
C．流程图：
开始 n=1; h=(b-a)/2; I1=0; I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h; 是 否 abs (I2–I1) >eps 否 n=n+1; h=(b-a) / n; I1=I2; I2=0; for i=0:n-1 是 x=a+h*i; x1=x+h; I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1)); I=I2; Step=n; 结束 D．应用举例
复合梯形法求数值积分应用举例，利用复合梯形法计算定积分a. 流程图
开始 421dx x21f= 1/(x^2-1) a=2; b=4 调用函数CombineTraprl(f,a,b,eps) 结束
b.程序代码：
[q,s]=CombineTraprl('1/(x^2-1)',2,4) %精度为默认的10-4
c.结果:
[q,s]=CombineTraprl('1/(x^2-1)',2,4,1.0e-5) %精度为10-5
结果：
所以从复合梯形公式可以得出
421dx0.2939 2x1二、1.将100个学生5门功课的成绩存入矩阵P中，进行如下处理： （1）分别求每门课的最高分、最低分及相应学生序号。 （2）分别求每门课的平均分和标准方差。
（3）5门课总分的最高分、最低分及相应学生序号。
（4）将5门课总分按从大到小顺序存入zcj中，相应学生序号存入xsxh。
A.流程图：
开始 构造一个矩阵做学生序号 再构造一个随机矩阵做学生成绩 让每门课最高分放入a矩阵，最低分放人b矩阵 使用函数max求出每门课最大值，并用find函数找到，然后输入到a； 同理用min求出每门课最小值。 使用mean函数求出每门课平均分c 使用var函数求出每门课的标准方差 使用sum函数求和，并找到每门课的最高分放入a2矩阵,每门课最低分放入b2。 用sort函数和descend把总分排序放人zcj文件和xsxh文件 结束 B.源代码
num=[1:1:100]';
p=[num,floor(50*rand(100,5)+45)]; a=[];b=[]; for i=1:5 m=max(p(:,i+1));
pp=find(p(:,i+1)==max(p(:,i+1))); aa=ones(size(pp))*m; a=[a;0,0;pp,aa] end for i=1:5 m=min(p(:,i+1));
pp=find(p(:,i+1)==min(p(:,i+1))); bb=ones(size(pp))*m; b=[b;0,0;pp,bb] end
c=mean(p(:,2:6)) d=var(p(:,2:6)) ps=sum((p(:,2:6))'); ps=ps';
a2=[find(ps==max(ps)),max(ps)] b2=[find(ps==min(ps)),min(ps)] [zcj,xsxh]=sort(ps,'descend')
C.结果：
（1）每门课的最高分、最低分及相应学生序号
（2）分别求每门课的平均分和标准方差。
（3）5门课总分的最高分、最低分及相应学生序号。
（4）将5门课总分按从大到小顺序存入zcj中，相应学生序号存入xsxh。
三、某气象观测站测得某日6：00~18：00之间每隔2h的室外温度（℃）如实验表1所示。
实验表1 室外温度观测结果（℃）
时间h 6 室内温18.0 度t1 室外温15.0 度t2 8 20.0 19.0 10 22.0 24.0 12 25.0 28.0 14 30.0 34.0 16 28.0 32.0 18 24.0 30.0 试用三次样条插值分别求该日室内外6：30～17：30之间每隔2h各点的近似温度（℃）
A．流程图
开始 x=[6 8 10 12 14 16 18]; y1=[18.0 20.0 22.0 25.0 30.0 28.0 24.0]; y2=[15.0 19.0 24.0 28.0 34.0 32.0 30.0]; x0=[6.30:2.0:17.30] 使用函数spline(x,y,x0)) 结束
B.代码：
x=[6 8 10 12 14 16 18]; %整体时间 y1=[18.0 20.0 22.0 25.0 30.0 28.0 24.0]; %室内温度 y2=[15.0 19.0 24.0 28.0 34.0 32.0 30.0]; %室外温度
x0=[6.30:2.0:17.30]; %外6：00～17：30之间每隔2h的时间 y11=spline(x,y1,x0) %使用三次样条插值函数 y12=spline(x,y2,x0) %使用三次样条插值函数
C．结果：
四、已知lgx在[1,101]区间10个整数采样点的函数值如下表所示，
试求lgx的5次拟合多项式p(x)，并绘制lgx和p(x)在[1,101]区间的函数曲线
A．流程图：
开始 x=[1 11 21 31 41 51 61 71 81 91 101]; lgx=[0 1.0414 1.3222 1.4914 1.6128 1.7076 1.7853 1.8513 1.9085 1.9590 2.0043]; 使用函数 polyfit (x,lgx,5) 结束
B.代码：
x=[1 11 21 31 41 51 61 71 81 91 101]; %10个整数采样点
lgx=[0 1.0414 1.3222 1.4914 1.6128 1.7076 1.7853 1.8513 1.9085 1.9590 2.0043]; %lgx 的函数值
[P,S] = polyfit(x,lgx,5); % lgx的5次拟合多项式
C.结果:
plot(x,lgx,':o') % X-lgx的曲线
2.521.510.50020406080100120
>> p(x) = polyval(P,x); % lgx的5次拟合多项式p(x) >> plot(x,p(x),'-*') % p(x)的函数曲线
2.521.510.50020406080100120
因篇幅问题不能全部显示，请点此查看更多更全内容