聪明文档网

西安交通大学计算方法B上机报告（完整题目和程序）

1. 计算以下和式：word/media/image1.gif word/media/image2.gif，要求：
(1)若保留11个有效数字，给出计算结果，并评价计算的算法；
(2)若要保留30个有效数字，则又将如何进行计算。

a) 解题思想

（1）先根据精度要求估计所需累加的项数n，使用后验误差估计方法，条件为：word/media/image3.gif (m为有效数字位数)。

（2）在该题中S的和式存在两个相近的数相减的问题，为了避免有效数字损失，在计算中改变了运算顺序，分别将正数和负数分别相加，然后再将其和相加。

（3）为了避免大数吃小数的问题，本题先计算出保留目标有效数字所需要的迭代次数，然后采用倒序相加的方法提高计算精度。

b) 算法实现的结构

1. S1=s2=0；

word/media/image4.gif;

2. for n=0,1,2,…,i

If word/media/image5.gif

end

3. for n=i,i-1,i-2,…,0

a1=4/(16^n*(8*n+1));

a2=2/(16^n*(8*n+4));

a3=1/(16^n*(8*n+5));

a4=1/(16^n*(8*n+6));

s1=a1+s1;

s2=a4+a3+a2+s2;

end

S=s1-s2;

c) 计算源程序

clear; %清除工作空间变量

clc; %清除命令窗口命令

m=input('请输入有效数字的位数m=');%输入需要的有效数字位数

S=0;s1=0;s2=0;%定义存储正数、负数和累加和的变量

for n=0:1:1000

t=(1/16^n)*(4/(8*n+1)-(2/(8*n+4)-1/(8*n+5)-1/(8*n+6)));

if t<=10^(-m) %根据有效数字位数确定所需累加的n值

break;

end

k=n;

end;

for n=(k-1):-1:0

a1=4/(16^n*(8*n+1));

a2=2/(16^n*(8*n+4));

a3=1/(16^n*(8*n+5));

a4=1/(16^n*(8*n+6));

s1=a1+s1; %第一项倒序累加

s2=a4+a3+a2+s2; %后三项倒序累加

end

S=s1-s2; %正数累加值和负数累加值的和

S=vpa(S,m) %控制S的精度

d) 计算结果与评价

当需保留11位有效数字时，需要将n值加到

n=7， s=3.1415926536；

当需保留30位有效数字时，需要将n值加到

n=22, s=3.14159265358979311599796346854。

由计算结果可以看出，采用从后往前进行计算的方式，避免了“大数吃小数”的问题，这种算法很好的保证了计算结果要求保留的准确数字有效位数的要求。另外，word/media/image6.gif中有多个负数相加，按照绝对值递增的顺序求和，减少了舍入误差带来的影响。

2. 某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：

(1)请用合适的曲线拟合所测数据点；

(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；

a) 解题思想

该题需要对所需光缆长度进行近似估计，通过对测量数据进行拟合并进行曲线积分即可得到河底光缆长度的近似值。由于数值点较多，如果使用多项式插值，会出现龙格现象导致误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但又希望将所得数据点都用上，且所用数据点越多越好，所以本题采用分段插值方式，即用分段多项式代替单个多项式作插值。分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略了水流对光缆的冲击，故本题使用三次样条插值进行计算精度是最高的。

b) 算法实现的结构

1 Forword/media/image7.gifword/media/image2.gif

1.1 word/media/image8.gif

2 Forword/media/image9.gif

2.1 For word/media/image10.gif

2.1.1 word/media/image11.gif

3 word/media/image12.gif

4 Forword/media/image13.gif

4.1 word/media/image14.gif

4.2 word/media/image15.gif

4.3 word/media/image16.gif

5 word/media/image17.gif

6 word/media/image18.gif

7 获取M的矩阵元素个数，记为m

8 Forword/media/image19.gif

8.1 word/media/image20.gif

8.2 word/media/image21.gif

8.3 word/media/image22.gif

9 word/media/image23.gif

10 For word/media/image24.gif

10.1 word/media/image25.gif

11 获取x的元素个数存入s

12 word/media/image26.gif

13 Forword/media/image27.gif

13.1 Ifword/media/image28.gifthenword/media/image29.gifbreak

elseword/media/image30.gif

14 word/media/image31.gifword/media/image32.gifword/media/image33.gif

word/media/image34.gif

15 计算光缆长度时，用如下公式：

word/media/image35.gif

word/media/image36.gif

c) 计算源程序

clear; %清除工作空间变量

clc; %清除命令窗口命令

x=0:1:20;%将分点位置以数组的形式存储于x中

X=0:0.2:20;%将插值后要求的点存储在X中

y=[9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.80,10.93];

%探测点位置的深度数据

n=length(x);%分点位置的个数

N=length(X);%利用插值函数待求位置的个数

%求三次样条插值函数s(x)，见课本P110

M=y;

for k=2:3;%计算三次样条的二阶差商并存放在数组M中

for i=n:-1:k;

M(i)=(M(i)-M(i-1))/(x(i)-x(i-k+1));

end

h(1)=x(2)-x(1);

%计算三对角阵系数a,b,c及右端向量d

for i=2:n-1;

h(i)=x(i+1)-x(i);

c(i)=h(i)/(h(i)+h(i-1));

a(i)=1-c(i);

b(i)=2;

d(i)=6*M(i+1);

end

M(1)=0;%补充自然边界条件

M(n)=0;

a(n)=0;

b(1)=2;

b(n)=2;

c(1)=0;d(1)=0;d(n)=0;

u(1)=b(1);%对三对角阵进行LU分解

Y(1)=d(1);

for k=2:n;

l(k)=a(k)/u(k-1);

u(k)=b(k)-l(k)*c(k-1);

Y(k)=d(k)-l(k)*Y(k-1);

end

M(n)=Y(n)/u(n);%TSS算法求解样条参数M(i)，见课本P49

for k=n-1:-1:1;

M(k)=(Y(k)-c(k)*M(k+1))/u(k);

end

s=zeros(1,N);

for m=1:N;

k=1;

for i=2:n-1

if X(m)<=x(i);

k=i-1;

break;

else

k=i;

end

H=x(k+1)-x(k);%在各区间用三次样条插值函数计算待求X点处的s值

x1=x(k+1)-X(m);

x2=X(m)-x(k);

s(m)=(M(k)*(x1^3)/6+M(k+1)*(x2^3)/6+(y(k)-(M(k)*(H^2)/6))*x1+(y(k+1)-(M(k+1)*(H^2)/6))*x2)/H;

end

%计算所需光缆长度

L=0;

for i=2:N

L=L+sqrt((X(i)-X(i-1))^2+(s(i)-s(i-1))^2);%利用第一类线积分求河底光缆的长度

end

disp('所需光缆长度为L=');

disp(L);

figure

plot(x,y,'*',X,s,'-')%绘制铺设河底光缆的曲线图

xlabel('分点位置','fontsize',16);%标注坐标轴含义

ylabel('测点深度/m','fontsize',16);

title('铺设河底光缆的插值曲线图','fontsize',16);

grid on;

d) 计算结果与评价

计算所得的运行结果：

铺设海底光缆的插值曲线图如所示：

word/media/image38.emf

上述插值结果表明，采用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势，可以有效避免龙格现象。在本题的计算中采用自然三次样条函数作为边界条件，在解线性方程组时使用了TSS算法求解带状矩阵，其计算速度快，是一种求解线性方程组的有效手段，在估计河底光缆长度时使用第一类线积分，最终计算出所需光缆的长度为 L=26.4844m。在计算编程过程中对三次样条插值的原理理解不透彻，导致编程逻辑思路不清晰，调试时间较长，体会到编程是一项仔细认真的工作，不能马虎！

3. 假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律;试计算这一天的平均气温,并试估计误差。

a) 解题思想

该题要求拟合一天的气温变化规律，由于记录数据点的数目较多。需要采用很高次的多项式做数据近似，这会给计算带来困难，也会导致误差很大。用“插值样条函数”做数据近似，虽然有很好的数值性质，且计算量也不大，但存放参数Mi的量很大，且没有一个统一的数学公式来表示，也带来了一些不便。另一方面，在有的实际问题中，用插值方法并不合适，当数据点的数目很大时，要求word/media/image39.gif通过所有数据点，可能会失去原数据所表示的规律，如果数据点是由测量而来的，必然带有误差，插值法要求准确通过这些不准确的数据点是不合适的。在这种情况下，不用插值标准而用选取word/media/image40.gif使word/media/image41.gif最小的标准，即采用最小二乘近似方法。

本题采用“最小二乘法”找出这一天的气温变化的规律，考虑使用最小二乘的二次函数、三次函数、四次函数以及指数函数word/media/image42.gif，计算相应的系数，估算误差，问题的难度是求解各种拟合函数的系数。由于利用法方程求解最小二乘系数时，方程的解不够稳定，所以本题利用正交化方法解最小二乘问题。

b) 算法结构

用正交化方法求数据的最小二乘近似问题，假定数据以用来生成了G，并将y作为其最后一列（第n+1列）存放，结果在a中，是误差是E2。

(一) 使用高次多项式拟合

1 将“时刻值”存入word/media/image43.gif，数据点的个数存入N

2 输入拟合多项式函数word/media/image44.gif的最高项次数n-1，则拟合多项式函数为word/media/image45.gif，。

3 根据给定数据点确定G,并将y作为n+1列存入G

Forword/media/image46.gif

Forword/media/image47.gif

word/media/image48.gifword/media/image49.gif

4 Forword/media/image50.gif

4.1 [形成矩阵word/media/image51.gif]

4.1.1 word/media/image52.gif

4.1.2 word/media/image53.gif

4.1.3 For word/media/image54.gif

4.1.3.1 word/media/image55.gif

4.1.4 word/media/image56.gif

4.2 [变换word/media/image57.gif到word/media/image58.gif]

4.2.1 word/media/image59.gif

Forword/media/image60.gif

4.2.2 word/media/image61.gif

4.2.3 Forword/media/image62.gif

4.2.3.1 word/media/image63.gif

5 [解三角方程word/media/image64.gif]

5.1 word/media/image65.gif

5.2 Forword/media/image66.gif

5.2.1 word/media/image67.gif

6 [计算误差word/media/image68.gif]

word/media/image69.gifword/media/image70.gif

(二) 使用指数函数拟合

当为指数函数时，假定指数函数的形式为word/media/image42.gif。只需将y值求对数，其主体部分不变，令word/media/image71.gif，则可得多项式: word/media/image72.gif。

c) MATLAB 源程序

clear; %清除工作空间变量

clc; %清除命令窗口命令

x=0:24; %存入时刻值

y=[15 14 14 14 14 15 16 18 20 20 23 25 28 31 34 31 29 27 25 24 22 20 18 17 16];

N=length(x);

n=input('请输入所需拟合函数的最高项次数n=');

plot(x,y,'r*')

grid;

hold on;

n=n+1;

G=zeros(N,n+1); %定义零矩阵

G(:,n+1)=y; %把y作为G的最后一列存放

C=zeros(1,n); %建立C矩阵来存放σ

E=0;

F=0; %定义计算用的中间变量

B=zeros(1,N); %建立B用来存放β

%根据已给数据生成矩阵G

for j=1:n

for i=1:N

G(i,j)=x(1,i)^(j-1);

end

%形成矩阵Qk，见课本P123

for k=1:n

for i=k:N

C(k)=G(i,k)^2+C(k);

end

C(k)=-sign(G(k,k))*(C(k)^0.5); %计算σ的值存入C中

w(k)=G(k,k)-C(k); %建立w来存放ω

for j=k+1:N

w(j)=G(j,k);

end

B(k)=C(k)*w(k); %计算β的值并存入B中

%变换矩阵Gk-1到Gk，参考课本P123

G(k,k)=C(k);

for j=k+1:n+1

E=0;

for i=k:N

E=w(i)*G(i,j)+E;

end

t=E/B(k);

for i=k:N

G(i,j)=t*w(i)+G(i,j);

end

%求解三角方程：Rx=h1

a(n)=G(n,n+1)/G(n,n);

for i=n-1:(-1):1

for j=i+1:n

F=G(i,j)*a(j)+F;

end

a(i)=(G(i,n+1)-F)/G(i,i); %a(i)存放各级系数

F=0;

end

disp(['各级系数分别为 ',num2str(a)])

%拟合后的回代求解过程

p=zeros(1,N);

for j=1:N

for i=1:n

p(j)=p(j)+a(i)*x(j)^(i-1);

end

plot(x,p,'bx',x,p,'-');

xlabel('时刻/h','fontsize',16);%标注坐标轴含义

ylabel('平均温度/℃','fontsize',16);

title('拟合平均气温曲线图','fontsize',16);

grid on;

E2=0; %用E2来存放误差

%计算拟合误差

for i=n+1:N

E2=G(i,n+1)^2+E2;

end

E2=E2^0.5;

disp(['拟合所得误差为E2=',num2str(E2)]);

Tm=0;

for i=1:N

Tm=Tm+p(i);

end

t=Tm/N;%计算平均温度

disp(['平均温度为t=',num2str(t),'℃'])

d) 计算结果与分析

（1）采用二次函数拟合

故可得函数形式为：word/media/image74.gif

二次函数的最小二乘曲线如下图所示：

word/media/image75.emf

（2）采用三次函数拟合

故可得函数形式为：word/media/image77.gif

三次函数的最小二乘曲线如下图所示：

word/media/image78.emf

（3）采用四次函数拟合

故可得函数形式为：word/media/image80.gif

四次函数的最小二乘曲线如下图所示：

word/media/image81.emf

（4）采用指数函数拟合

故可得函数形式为：word/media/image83.gif

指数函数的最小二乘曲线如下图所示：

word/media/image84.emf

通过上述拟合结果可以发现，使用最小二乘法进行数据拟合时，多项式的次数越高，拟合曲线也更加接近数据点，计算拟合的效果越好，误差越小，结果越准确。同时，指数多项式拟合的次数虽然不高，但误差最小，是因为指数型二次函数是对数据点求对数进行拟合的，故结果最准确。在进行最小二乘法编程时发现，使用正交化方法求解最小二乘问题比使用法方程方法更稳定，结果更可靠。但编程是中间变量较多容易造成混乱，计算量也较大，但本题在编程过程中有一定的适用性，只需要对数据点进行处理就可以拟合出给定次数的多项式，为以后的数据拟合工作打下了基础。

4. 设计算法，求出非线性方程word/media/image85.gif的所有实根，并使误差不超过word/media/image86.gif。

a) 解题思想

由题目所给的待求方程可知该方程为5次多项式，导数的求值比较容易，可采用牛顿法迭代求解。具体实现步骤为：首先，研究函数的形态，确定根的范围；通过剖分区间的方法确定根的位置，然后利用牛顿法的基本原理进行求解，找到满足精度要求的解。

令word/media/image87.gif，其迭代格式为word/media/image88.gifword/media/image2.gif

要使牛顿迭代法收敛，则必须寻找根的合理区间[a,b],使得在该区间内word/media/image89.gif，即有根。在选定的区间内函数word/media/image90.gif的一二阶导数word/media/image91.gif均不变号。选定一个初值word/media/image92.gif，使word/media/image93.gif，则牛顿迭代求解过程完成。

b) 算法结构

1 Forword/media/image94.gif

1.1 word/media/image95.gif

1.2 Ifword/media/image96.gifthen

1.2.1 输出“word/media/image90.gif的值已足够小”

1.2.2 word/media/image97.gif；stop

else

1.2.3 word/media/image98.gif

1.2.4 word/media/image99.gif

1.2.5 Ifword/media/image100.gifthen

1.2.5.1 输出“word/media/image101.gif已足够小”；stop

else

1.2.5.2 word/media/image102.gif

2 输出“word/media/image103.gif次迭代后仍不收敛”，stop

c) 计算源程序

clear;

clc;

%根据函数曲线观察根的大致分布区间

X=-8:0.1:8;

for i=1:161

Y(i)=6*X(i)^5-45*X(i)^2+20;

end

plot(X,Y,'r-','LineWidth',2);%画出Y随X变化的曲线图

%标注坐标轴含义及图像标题

xlabel('X','fontsize',10);

ylabel('Y','fontsize',10);

title('多项式函数变化曲线图','fontsize',10);

grid on;

k=1;

Y(162)=0;

%使用二分法来确定不同根区间的迭代初值，并存于a中

for i=1:161

if Y(i)*Y(i+1)<0

a(k)=X(i);

k=k+1;

continue

end

x=a;

n=length(x);

%%%牛顿迭代法计算方程根?

E=1e-4;

for i=1:n

for j=1:50

f=x(i)-(6*x(i)^5-45*x(i)^2+20)/(30*x(i)^4-90*x(i));

b=abs(f-x(i));

if b判断误差是否满足要求

break

end

x(i)=f;%求得的根重新赋值给x

end

disp(['计算结果为:',num2str(x)]);

d) 计算结果与分析

（1）根的大致分布区间曲线图如下

word/media/image104.emf

（2）根的求解结果：

（3）分析总结

由计算结果可知三个根分别为：-0.65455, 0.68118, 1.870，满足误差精度要求，在计算求解该题过程中，利用已知多项式函数的图像可以初步观察出根的大致分布区间，然后通过循环得出根区间的迭代初值，从而有效缩减了求根区间。为防止迭代进入死循环，设置了最高迭代次数为50次，而结果显示方程的三个根的迭代次数均在4次以内，可以看出牛顿迭代法的迭代求解多项式根的速度非常快，精度也高。

5. 线性方程组求解。

（1）编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。所附方程组的类型为对角占优的带状方程组。

（2）针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。

(一)

a) 解题思想

高斯消去法的思想很简单，它的第一步是将n元方程组的n-1个方程实施“消元”，形成一个与原方程等价的新方程组，而被消元的n-1个方程实际上形成了一个n-1元方程组，逐步消元后得到一个上三角形状的方程组，它的从上到下的各个方程逐个减少一个未知量，最后一个方程是一元一次方程，然后，从最后一个方程解出一个未知量开始，逐次回代求得方程组的全部解。

列主元消去法是当高斯消元到第k步时，从k列的word/media/image106.gif以下（包括word/media/image106.gif）的各元素中选出绝对值最大的，然后通过行交换将其交换到word/media/image106.gifword/media/image2.gif的位置上，这样就可以有效的避免大数除以小数而造成的上溢现象。根据教材提供的算法，编写高斯列主元消去法程序，编写列主元消去法的子函数与适应于超大规模超出系统内存的方程组的改编程序。

b) 算法结构

1.数据文件的文件名后缀为.dat,形式为：文件名.dat；

2.数据文件中的数据均为二进制记录结构，因此必须使用二进制方式进行读取；

3.数据文件的结构，分为以下四个部分：

(1)文件标示部分,该部分用于存放数据文件的描述信息

结构如下(用C语言格式进行描述)：

typedefstructFileInfo {

long int id; // 数据文件标示

long intver; // 数据文件版本号

long int id1; // 备用标志

} FILEINFO;

其中：

①id:为该数据文件的标识，值为0xF1E1D1A0即为：十六进制的F1E1D1A0

②ver:为数据文件的版本号，值为16进制数据,

③id1:为备用标志字段，暂时未用;

(2)矩阵描述部分:此部分中包括矩阵的阶数和上下带宽,如果是稀疏矩阵，则上下带宽值为0。

结构如下：

typedefstructHeadInfo {

long int n; // 方程组的阶数

long int q; // 带状矩阵上带宽

long int p; // 带状矩阵下带宽

} HEADINFO;

(3) 系数矩阵数据部分：该部分存放方程组系数矩阵中的所有元素

①如存贮格式为非压缩格式,则按行方式顺序存贮系数矩阵中的每一个
元素，元素总个数为n*n,每个元素的类型为float型；

②如果存贮格式是压缩方式,则同样是按行方式进行存贮,每行中只放上下带宽内的非零元素,即每行中存贮的元素都为p+q+1个；

(4)右端系数部分：该部分存放方程组中的右端系数

按顺序存贮右端系数的每个元素，个数为n个，每个系数的类型为float型

4 数据文件的读取

f=fopen('fun003.dat','r'); %打开文件，.dat文件放在 m 文件同一目录下，

a=fread(f,3,'uint') %读取头文件，3-读取前 3 个，若读取压缩格式的，头文件为 5 个

b=fread(f,inf,'float'); %读取剩下的文件，float 型

id=dec2hex(a(1));

ver=dec2hex(a(2)); %这两句是进行进制转换，读取 id 与ver

5 A的阶数word/media/image107.gif

6 Forword/media/image108.gif

6.1 找满足word/media/image109.gif的下标word/media/image110.gif

6.2 Forword/media/image111.gif

6.2.1 word/media/image112.gif

6.3 word/media/image113.gif

6.4 Forword/media/image114.gif

6.4.1 word/media/image115.gif

6.4.2 Forword/media/image116.gif

6.4.2.1 word/media/image117.gif

6.4.3 word/media/image118.gif

word/media/image119.gif

For word/media/image120.gif

word/media/image121.gif

c) 计算源程序

clear; %清除工作空间变量

clc; %清除命令窗口命令

%%读取系数矩阵

[f,p]=uigetfile('*.dat','选择数据文件'); %读取数据文件

num=5; %输入系数矩阵文件头的个数

name=strcat(p,f);

file=fopen(name,'r');

head=fread(file,num,'uint'); %读取二进制头文件

id=dec2hex(head(1)); %读取标识符

fprintf('文件标识符为');

ver=dec2hex(head(2)); %读取版本号

fprintf('文件版本号为');

ver

n=head(3); %读取阶数

fprintf('读取矩阵的阶数');

q=head(4); %上带宽

fprintf('读取矩阵的上带宽');

p=head(5); %下带宽

fprintf('读取矩阵的下带宽');

dist=4*num;

fseek(file,dist,'bof'); %把句柄值转向第六个元素开头处

[A,count]=fread(file,inf,'float'); %读取二进制文件，获取系数矩阵

fclose(file); %关闭二进制头文件

%%对非压缩带状矩阵进行求解

ifver=='102'

a=zeros(n,n);

for i=1:n

for j=1:n

a(i,j)=A((i-1)*n+j); %求系数矩阵a(i,j)

end

b=zeros(n,1);

for i=1:n

b(i)=A(n*n+i);

end

for k=1:n-1 %列主元高斯消去法

m=k;

for i=k+1:n, %寻找主元

if abs(a(m,k))

m=i;

end

if a(m,k)==0 %遇到条件终止

disp('there is a mistake！')

return

end

for j=1:n %交换元素位置的主元

t=a(k,j);

a(k,j)=a(m,j);

a(m,j)=t;

t=b(k);

b(k)=b(m);

b(m)=t;

end

for i=k+1:n %计算l(i,k)并将其放到a(i,k)中

a(i,k)=a(i,k)/a(k,k);

for j=k+1:n

a(i,j)=a(i,j)-a(i,k)*a(k,j);

end

b(i)=b(i)-a(i,k)*b(k);

end

x=zeros(n,1); %回代过程

x(n)=b(n)/a(n,n);

for k=n-1:-1:1

x(k)=(b(k)-sum(a(k,k+1:n)*x(k+1:n)))/a(k,k);

end

%%对压缩带状矩阵进行求解

if ver=='202' %高斯消去法

m=p+q+1;

a=zeros(n,m);

for i=1:1:n

for j=1:1:m

a(i,j)=A((i-1)*m+j); %求a(i,j)

end

b=zeros(n,1);

for i=1:1:n

b(i)=A(n*m+i); %求b(i)

end

for k=1:1:(n-1) %开始消去过程

if a(k,(p+1))==0

disp('there is a mistake！');

break;

end

st1=n;

if (k+p)

st1=k+p;

end

for i=(k+1):1:st1

a(i,(k+p-i+1))=a(i,(k+p-i+1))/a(k,(p+1));

for j=(k+1):1:(k+q)

a(i,j+p-i+1)=a(i,j+p-i+1)-a(i,k+p-i+1)*a(k,j+p-k+1);

end

b(i)=b(i)-a(i,k+p-i+1)*b(k);

end

x=zeros(n,1); %回代过程

x(n)=b(n)/a(n,p+1);

sum=0;

for k=(n-1):-1:1

sum=b(k);

st2=n;

if (k+q)

st2=k+q;

end

for j=(k+1):1:st2

sum=sum-a(k,j+p-k+1)*x(j);

end

x(k)=sum/a(k,p+1);

sum=0;

end

disp('方程组的的解为：') %输出方程组的解

disp(x)

d) 计算结果与分析

（1） dat51.dat求解结果：

文件标识符为id =F1E1D1A0

文件版本号为ver =102

读取矩阵的阶数n =15

读取矩阵的上带宽q =3

读取矩阵的下带宽p =3

方程组的的解为：

1.0000

（2） dat52求解结果

文件标识符为id =F1E1D1A0

文件版本号为ver =202

读取矩阵的阶数n =20

读取矩阵的上带宽q =5

读取矩阵的下带宽p =5

方程组的的解为：

1.0000

（3） dat53求解结果

文件标识符为id =F1E1D1A0

文件版本号为ver =102

读取矩阵的阶数n =2160

读取矩阵的上带宽q =5

读取矩阵的下带宽p =5

方程组的的解为：

3.1400

word/media/image122.gif

3.1399

3.1402

3.1400

3.1404

3.1396

3.1402

3.1392

3.1375

3.1406

3.1400

3.1403

3.1483

3.1487

3.1689

3.1529

3.1910

3.0728

3.0698

2.9621

3.3051

2.9243

3.5790

3.8189

3.5284

4.1508

4.7523

-1.6718

1.3798

4.9003

16.6720

4.2458

2.1043

-4.1581

0.9427

6.8407

2.3002

4.5557

2.5026

2.7366

3.0056

3.3210

2.9951

3.1727

3.1951

3.1660

3.1196

3.1558

3.1567

3.1386

3.1396

3.1420

3.1302

3.1379

3.1401

3.1418

3.1416

3.1404

3.1401

3.1397

3.1394

3.1402

3.1401

3.1400

word/media/image122.gif

3.1400

（4） dat54求解结果

文件标识符为id =F1E1D1A0

文件版本号为ver =202

读取矩阵的阶数n =43240

读取矩阵的上带宽q =4

读取矩阵的下带宽p = 4

方程组的的解为：

5.2100

word/media/image122.gif

5.2100

(二) 实际问题算例（摘自数值传热学作业）

1. 一正方形导热物体的底边是绝热的，其余三边上的温度如下图所示，试确定该正方形内节点1、2、3、4的温度，导热物体无内热源，物性为常数。

解：导热物体无内热源，物性为常数的控制方程可以写成

word/media/image124.gif

取均分网格，word/media/image125.gif，采用中心差分进行离散可得如下离散方程：

word/media/image126.gif

因为无内热源，均分网格，所以可取word/media/image127.gif，则word/media/image128.gif。

即有：word/media/image129.gif故温度待求节点的方程整理为：

word/media/image130.gif

word/media/image131.gif

word/media/image132.gif

word/media/image133.gif

word/media/image135.gif-4 1 1 0 word/media/image137.gif-70word/media/image134.gifword/media/image136.gif

1 -4 0 1 word/media/image138.gifword/media/image139.gif -50

1 0 -3 1 word/media/image140.gif-15

0 1 1 -3 word/media/image141.gif-10

用列主元高斯消去法对上述矩阵进行求解，其程序如下：

clear; %清除工作空间变量

clc; %清除命令窗口命令

A=[-4 1 1 0 -70;

1 -4 0 1 -50;

1 0 -3 1 -15;

0 1 1 -3 -10];

[m,n]=size(A); %计算矩阵的维数

%将绝对值最大系数作为主元

fori=1:(m-1)

max=i;

fork=i+1:m

ifabs(A(k,i))>abs(A(max,i))

max=k;

end

temp=A(i,i:m);

A(i,i:m)=A(max,i:m);

A(max,i:m)=temp;

for j=(i+1):m

A(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i);

end

%消元完成后进行回代求解x

disp('进行回代求解')

x(m)=A(m,n)/A(m,m);

for k=(m-1):-1:1

x(k)=(A(k,n)-A(k,k+1:m)*x(k+1:m)')/A(k,k);

end

disp(['求解的结果为：',num2str(x)])

2. 计算结果与分析

故T1=28.7368，T2=24.2632，T3=20.6842，T4=18.3158。运算结果和正确结果吻合。

对于附件中所给的方程组，求解的根都相等，其中dat51.dat，为15阶方程，根均为1.0000；dat52.dat为20阶，其根均为1.0000；dat53.dat为2160阶，其根不均为3.1400, 具体结果上述已提及；dat54.dat为43240阶，其根均为5.2100。在求解非压缩格式的方程根时，采用一般列主元高斯消去法，计算量大，耗时长。在求解压缩格式的方程根时，由于方程阶数较高，运算时间太长，对程序进行了改进，只对非零元素进行计算，运算步数大量减少，速度得到了极大的提高，采用此方法求解方程组非常高效便捷，尤其在求解高阶方程组过程中其优越性更加明显。

感悟：

非常开心选了老师的《计算方法》这门课程，老师课堂上详细认真的粉笔板书让我印象深刻，也让我对数值计算产生了极大的兴趣，对此门课程有了一个全新的识。

老师您给我们布置的上机作业，让我将理论上的东西运用到了实际计算应用中，也是对自己知识掌握程度的一次极好的检验，从开始的算法思路，到编写源程序，再到运行调试，以及最后运行出正确结果，对于我来说是一个很好的学习和锻炼的过程。不仅对我学到的理论知识进行了巩固，也培养了分析解决实际问题的能力。使我感受到应用所学知识解决问题的成绩感和喜悦感，真的是受益匪浅，最后谢谢老师知识的传授！。

《西安交通大学计算方法B上机报告（完整题目和程序）.doc》

将本文的Word文档下载，方便收藏和打印