计算方法B上机报告

西安交通大学

计算方法B上机报告

班级：XXXXXXXX 姓名：XXX

学号：XXXXXXXX

2015/12/13 Sunday

目录

目录

1 题目一 ............................................................................................................................... 1 1. 数值计算........................................................................................................................ 1 2. 实现思想........................................................................................................................ 1 3. 源程序............................................................................................................................ 1 4. 计算结果........................................................................................................................ 2 5. 分析总结........................................................................................................................ 2

2 题目二 ............................................................................................................................... 3 1. 数据近似........................................................................................................................ 3 2. 实现思想........................................................................................................................ 3 3. 源程序............................................................................................................................ 5 4. 计算结果........................................................................................................................ 6 5. 分析总结........................................................................................................................ 7

3 题目三 ............................................................................................................................... 8 1. 数据拟合........................................................................................................................ 8 2. 实现思想........................................................................................................................ 8 3. 源程序............................................................................................................................ 9 4. 计算结果...................................................................................................................... 12 5. 分析总结...................................................................................................................... 14

4 题目四 ............................................................................................................................. 15 1. 非线性方程求解 .......................................................................................................... 15 2. 实现思想...................................................................................................................... 15 3. 源程序.......................................................................................................................... 16 4. 计算结果...................................................................................................................... 17 5. 分析总结...................................................................................................................... 17

5 题目五 ............................................................................................................................. 18 1. 线性方程组求解。 ...................................................................................................... 18 2. 文件格式说明 .............................................................................................................. 18 3. 实现思想...................................................................................................................... 19 4. 源程序.......................................................................................................................... 20

1

目录

5. 计算结果...................................................................................................................... 23 6. 分析总结...................................................................................................................... 25 7. 实际问题...................................................................................................................... 25

2

计算方法B上机报告――题目一

1 题目一

1. 数值计算

计算以下和式：S??1nn?016?211??4????8n?18n?48n?58n?6?，要求：

??(1)若保留11个有效数字，给出计算结果，并评价计算的算法；

(2)若要保留30个有效数字，则又将如何进行计算。

2. 实现思想

对于(1)可用迭代的方法进行处理。通过do-while循环使和式从0开始计算直到结果满足有效数字即可。在循环中通过比较本次和上一次结果之差的绝对值与相应有效数字的大小的1/2，即可构成循环中止条件。

对于(2)可用同样的方法进行计算，然而由于所保留的有效数字较大，此时若不对上述算法进行改善的话，由于每一步所得计算结果都是其真实值，此时可能会存在有效数字丢失的情况，从而影响计算的准确性。因此我们需要对程序中的计算精度进行控制。在Maltlab中可利用digits()函数对运算精度进行规定，为防止出现有效位的损失，在实际程序中需将精度提高几位。而对迭代中每一步所得到的结果则利用vpa()函数进行限定，这样一来我们对每一个运算都控制了精度，而不只是单纯控制了结果的精度

3. 源程序

使用Matlab所编写程序如下所示：

(1)

clear;clc;

sd = 11; %所要保留有效数字 n = 0; %迭代次数

S1 = 0; %当前n时计算结果 S2 = 0; %n-1时计算结果 while 1

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

if abs(S1-S2)<0.5*10^(-(sd-1)) %迭代是否终止判断条件 break end

S2 = S1; n = n+1; end

1

计算方法B上机报告――题目一

S=vpa(S1,sd) (2)

clear;clc;

sd = 30; %所要保留有效数字 digits(sd+2);%控制精度 n = 0; %迭代次数

S1 = vpa(0); %当前n时计算结果 S2 = vpa(0); %n-1时计算结果 while 1

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

if abs(S1-S2)

S2 = S1; n = n+1; end

S=vpa(S1,sd)

4. 计算结果

(1)

S =3.1415926536 (2)

S =3.14159265358979323843571167922

5. 分析总结

从计算结果可以看出与准确值3.1415926535897932384357116792180407…相比(1)和(2)所得结果分别是准确值保留11位和30位有效数字的结果，即利用本程序所得结果准确。

而如果直接利用程序(1)对(2)进行求解，即只是单独修改有效数字位数sd而没有对运算精度加以限制所得结果为：3.14159265358979323846264338328。我们可以看出该结果与程序(2)所得结果以及准确值相差较大。这是由于未对计算过程中的每一步运算的运算精度进行控制造成有效数字的损失。当所需保留的有效数字较低时，该影响尚不明显，而当计算量较大时会对计算结果产生较大影响，使最终结果和准确值偏离较大。因此，在进行运算所要求保留的有效数字较大时要对计算结果的运算精度进行控制，防止运算结果产生较大误差。

2

计算方法B上机报告――题目二

2 题目二

1. 数据近似

某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示： 分点 深度 分点 深度 分点 深度 0 9.01 7 1 8.96 8 2 7.96 9 3 7.97 10 4 8.02 11 5 6 9.05 10.13 12 13 11.18 12.26 13.28 13.32 12.61 11.29 10.22 14 9.15 15 7.90 16 7.95 17 8.86 18 19 20 9.81 10.80 10.93 (1)请用合适的曲线拟合所测数据点；

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

2. 实现思想

本题首先要对所测量的21个点进行数据拟合，再对拟合的曲线进行线积分运算即可得到所需光缆长度的近似值。

数据近似的方法一般有两种，一个是多项式插值，另一种是分段插值。如果将所有的测量点都用上，采用多项式插值时，由于龙格现象的存在，随着点数的增大，有时会在两端产生激烈的震荡，使函数不收敛。因此，为了将所有的测量点都用上，本题应采用分段插值方法求解。而为了改善分段插值所引起的曲线的光滑性不够的问题，常采用较高的分段多项式改善曲线的光滑性质。故本题采用三次样条插值。

三次样条插值算法结构如下所示： 算法SPLINEM(x,y,M,n,?0,d0,?n,dn)

1.

For i?0,1,2,?,n

1.1 yi?Mi

2. For k?1,2

3

计算方法B上机报告――题目二

2.1 For i?n,n?1,?,k

2.1.1

3. 4.

(Mi?Mi?1)/(xi?xi?k)?Mi

xi?x0?h1

For i?0,1,2,?,n?1

4.1 xi?1?xi?hi?1

4.2 hi?1/(hi?hi?1)?ci;1?ci?ai;2?b 4.3 6Mi?di

5.

d0?M0;dn?Mn;?0?c0;2?b0;?n?an;2?bn

6. 将M矩阵中的元素个数存入m中 7. 8.

b1??1,d1??1

For k?2,3,?,m

8.1 ak/?k?1?lk 8.2 bk?lk?ck?1??k 8.3 dk?lk??k?1??k

9. ?m/?m?Mm

10. For k?m?1,m?2,?,1

10.1 (?k?ck?Mk+1)/?k?Mk 11. 将x中的元素个数存入r中 12. 1?k

13. For i?0,1,2，?,r?1

??xi then i?k; break 13.1 if xelse i?1?k;

??x;x??xk?1?x 14. xk?xk?1?h;xk?x?4

计算方法B上机报告――题目二

[Mk-1?x3x3h2h2?? ?Mk?(yk?1?Mk?1)x?(yk?Mk)x]/h?y6666其中Mi记在xi处的二阶导数s''(xi)，在自然样条插值中?0?0,?n?0,d0?0,dn?0。

得到河底光缆曲线的函数s(x)后，可利用第一类线积分计算光缆长度l为：

l??ds??0202001?(s(x))dx???'2k?019k?1k1?(s'(x))2dx

3. 源程序

按照上述算法所编写的MATLAB程序如下所示：

clear;clc;

x=0:1:20; %等分点数组 X=0:0.2:20; %插值点数组

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)

M=y; %计算y的二阶差商并存储在矩阵M中 for k=2:3;

for i=n:-1:k;

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

%计算三对角阵系数a,b,c及右端向量d h(1)=x(2)-x(1); 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;b(1)=2;b(n)=2;c(1)=0;a(n)=0; d(1)=0;d(n)=0; %追赶法计算样条参数M u(1)=b(1); ynew(1)=d(1); %追

for k=2:n;

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

u(k)=b(k)-len(k)*c(k-1); ynew(k)=d(k)-len(k)*ynew(k-1);

5

计算方法B上机报告――题目二

end

M(n)=ynew(n)/u(n); %赶

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

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

for m=1:N; k=1;

for i=2:n-1

if X(m)<=x(i); k=i-1; break; else k=i; end end

%在各区间用三次样条插值函数计算X点处的值 h=x(k+1)-x(k);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

%计算所需光缆长度

Length=0; for i=2:N

Length=Length+sqrt((X(i)-X(i-1))^2+(s(i)-s(i-1))^2); end

%河底光缆长度 Length

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

xlabel('宽度');ylabel('深度/m');title('铺设河底光缆的曲线图');grid;

4. 计算结果

采用三次样条插值拟合后的图形如下所示：

6

计算方法B上机报告――题目三

%计算误差e

rou=sum(G(n+1:m,n+1).^2); %结果显示

fprintf('\\n指数函数的系数是：a=%f,b=%f,c=%f',a,b,c); fprintf('\\n使用指数函数拟合的误差是:%f\\n',rou); pexp=a.*exp(-b.*(x-c).^2); for i=1:m

pexp(i)=exp(a+b*x(i)+c*x(i)^2); end

figure(2);

%绘制指数函数拟合曲线 plot(x,y,'k.-',x,pexp,'r*-'); title(['指数函数拟合曲线']);

legend('实测平均温度温度曲线','指数函数拟合曲线',n) xlabel('时刻/h');ylabel('平均气温/℃');xlim([0,24]); grid;

%计算平均温度 Tav=sum(y(1:m))/m;

disp(['平均温度为',num2str(Tav),'℃']);

4. 计算结果

(1) 二次多项式函数拟合程序运行结果如下所示：

2次多项式函数的系数由低到高为: 8.3063 2.6064 -0.0938 使用2次多项式函数拟合的误差是280.339547: 平均温度为21.2℃

(2) 三次多项式函数拟合程序运行结果如下所示：

12

计算方法B上机报告――题目三

3次多项式函数的系数由低到高为: 13.3880 -0.2273 0.2075 -0.0084

使用3次多项式函数拟合的误差是131.061822: 平均温度为21.2℃

(3) 四次多项式函数拟合程序运行结果如下所示：

4次多项式函数的系数由低到高为: 16.7939 -3.7050 0.8909 -0.0532 0.0009

使用4次多项式函数拟合的误差是59.041189: 平均温度为21.2℃

13

计算方法B上机报告――题目三

(4) 指数函数p(x)?e?0??1x??2x拟合程序运行结果如下所示：

指数函数的系数是：a=2.383470,b=0.125093,c=-0.004442 使用指数函数拟合的误差是:57.034644 平均温度为21.2℃

2

5. 分析总结

通过比较不同多项式的拟合效果可以看出，随着多项式次数的增加曲线拟合效果越来越好，误差也逐渐减小。因此，为了得到更好的拟合效果，应尽量采用较高次多项式。当数据具有指数的性质时，采用指数函数更为合适，可在比多项式阶数更少的情况下获得更小的误差。

14

计算方法B上机报告――题目四

4 题目四

1. 非线性方程求解

设计算法，求出非线性方程6x5?45x2?20?0的所有实根，并使误差不超过

10?4

2. 实现思想

本题要对非线性方程进行求解，可采用迭代法进行计算。迭代法包括简单迭代法、牛顿法、割线法以及区间法。虽然牛顿迭代法、割线法的收敛速度较快，但容易引起收敛性问题，以及有时方程计算结果虽然满足给定误差界但解和准确解的误差依然有可能很大的情况。因此我们这里采用区间法进行求解。因为区间法总是收敛的，且它的迭代过程确定了一个区间序列，使得每个区间都包含方程的某个解，只要区间长度小于给定误差界时，解和准确解的误差就满足所需误差要求。

常用的区间方法为二分法，每次迭代使区间长度减小一半，其算法结构如下所示：

算法BISECTION(x(0),x(0),f,?1,?2,x)

1. 2. 3. 4.

f(x(0))?f0;f(x(1))?f1

Iff0f1?0then stop

Iff0??2then输出x(0)作为根; stop Iff1??2then输出x(1)作为根; stop 1(0)x?x(1)??x ?2Ifx(1)?x??1x(1)then输出x作为根; stop

5. 6. 7. 8.

f(x)?f

Iff??2then输出x作为根;

15

计算方法B上机报告――题目四

9.

Iff0f1?0then

9.1 x?x(0);f?f0 else

9.2 x?x(1);f?f1

10. go to 5

算法中?1表示相邻两侧迭代的差的误差范围，?2表示某次迭代的方程的计算结果的误差范围。依题意本题只考虑相邻两侧迭代的差的误差，即解的误差。

3. 源程序

按照上述算法所编写的MATLAB程序如下所示：

clc;clear;

%确定根所在的区间 interval=[];

for i=-10:0.1:10

if (fx(i)*fx(i+0.1)<0) interval=[interval i]; end end

%确定根的个数 n=length(interval);

%对每个根区间采用二分法求解 for j=1:n

x0=interval(j); %左端点 x1=interval(j)+1;%右端点

while (x1-x0)>10^(-4) if fx(x0)*fx((x0+x1)/2)>0 x0=(x0+x1)/2; else

x1=(x0+x1)/2; end end

root(j)=x0; end

%方程组误差不超过10^-4的解 root

其中程序中的fx函数为：

function [y]=fx(x) y=6*x^5-45*x^2+20; end

16

计算方法B上机报告――题目四

4. 计算结果

程序运行结果为：

root= -0.6546 0.6811 1.8707

对结果进行验证，分别将-0.6546与-0.6545、0.6811与0.6812、1.8707与1.8708代入fx函数中得。

fx(-0.6546)= -0.0031 fx(0.6811)= 0.0032 fx(1.8707)= -0.0118

fx(-0.6545)= 0.0034 fx(0.6812)= -0.0023 fx(-1.8708)= 0.0081

由计算结果可以看出所得解在相差10^-4时的函数的结果均为一正一负，即所得解满足题目所要求的解的误差范围。

5. 分析总结

利用区间方法求解非线性方程的解时虽然其收敛速度不如牛顿迭代和割线法，但其保证了算法的收敛性，保证了解和准确解之差满足误差范围。

但从计算结果中我们同样可以看出利用区间法得到的解回代到原方程中所得解与0最大的误差相差0.01以上，如果本题还要求利用所求解运算后的方程在一定误差范围内时，此时虽然在判断条件上加上计算结果的误差限定同样得到误差范围内的解，但此时利用区间法的计算量将很大，不如采用牛顿法或割线发经济。因此，处理不同求解非线性方程根的问题时，应具体问题具体分析，兼顾收敛速度与收敛性，选择合适的迭代方法。

17

计算方法B上机报告――题目五

5 题目五

1. 线性方程组求解。

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

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

2. 文件格式说明

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

(2) 数据文件中的数据均为二进制记录结构，因此必须使用二进制方式进行读取； (3) 数据文件的结构，分为以下四个部分：

a) 文件标示部分,该部分用于存放数据文件的描述信息结构如下(用C语言格

式进行描述)：

typedefstructFileInfo { long int id; // 数据文件标示 long intver; // 数据文件版本号 long int id1; // 备用标志

} FILEINFO; 其中：

① id:为该数据文件的标识，值为0xF1E1D1A0即为：十六进制的F1E1D1A0 ② ver:为数据文件的版本号，值为16进制数据, 版本号 0x101 0x102 0x201 0x202

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

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

上下带宽值为0。 结构如下：

typedefstructHeadInfo {

long int n; // 方程组的阶数 long int q; // 带状矩阵上带宽

说明 系数矩阵为非压缩格式稀疏矩阵 系数矩阵为非压缩格式带状矩阵 系数矩阵为压缩格式稀疏矩阵 系数矩阵为压缩格式带状矩阵 18

计算方法B上机报告――题目五

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

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

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

② 如果存贮格式是压缩方式,则同样是按行方式进行存贮,每行中只 放上下带宽内的非零元素,即每行中存贮的元素都为p+q+1个； d) 右端系数部分：该部分存放方程组中的右端系数

按顺序存贮右端系数的每个元素，个数为n个，每个系数的类型为float型 (4) 数据文件说明：

a) Dat51.dat 为非压缩带状方程组，阶数为15阶，该方程组供调试程序

使用,该方程组的根都为1;

b) Dat52.dat 为压缩带状方程组，阶数为20阶，该方程组供调试程序使

用,该方程组的根都为1;

c) Dat53.dat 为非压缩带状方程组，阶数在2000阶左右; d) Dat54.dat 为压缩带状方程组，阶数在40000阶左右;

3. 实现思想

首先,通过利用MATLAB中自带的文件操作函数对dat文件进行处理，从而获取文件信息。对压缩带状方程组和非压缩带状方程组分别采用不同方法进行处理。

消去法的中心是降维，即将求解的n元方程组的问题转化为先解n-1元方程组，一旦次n-1元方程组的解取得，则剩余的一个未知量自然可以求得。通过这样逐步减少未知量的个数，从而对多元方程组进行求解。消去法包含两个主要步骤：消去和回代。在消去的过程中通过适当的交换方程的顺序保证消去过程能够顺利进行及计算解的精确度，交换的原则是使第k步消去过程中产生的数中绝对值最大的一个换到(k,k)位置而成为第k步消去的主元,即列主元Guss消去法。 其算法结构如下所示： 算法GAUSSPP(A,b,n,x)

1.

For k?1,2,3,?,n?1

1.1 找满足??kk=max?ik的下标?

i?kk1.2 If??kk?0then输出错误信息;stop 1.3 For j?1,2,3,?,n

1.3.1

?kj???j

k1.4 ?k???k

19

计算方法B上机报告――题目五

1.5 For i?k?1,k+2,?,n

1.5.1 1.5.2

?ik?kk??ik

For j?k?1,k+2,?,n

1.5.2.1 ?ij??ik?kj??ij 1.5.3

2. [回代过程]

2.1

?i??ik?k??i

?n?nn?xn

2.2 For k?n?1,n?2,?,1

2.2.1

n?????x?k?kjj?/?kk?xk

j?k?1??4. 源程序

按照上述算法所编写的MATLAB程序如下所示：

clc;clear;

%文件操作%读取系数矩阵

[filename,pathname]=uigetfile('*.dat','选择数据文件');%选择读取数据文件 num=5;%输入系数矩阵文件头的个数

name=strcat(pathname,filename);%连接字符构造文件路径 file=fopen(name,'r');

head=fread(file,num,'uint');%读取二进制头文件 id=dec2hex(head(1));%读取标识符 disp('文件标识符为');id

ver=dec2hex(head(2));%读取版本号 disp('文件版本号为');ver n=head(3); %读取阶数 disp('矩阵A的阶数');n q=head(4); %上带宽

disp('矩阵A的上带宽');q p=head(5);%下带宽

disp('矩阵A的下带宽');p dist=4*num;

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

[A,count]=fread(file,inf,'float');%读取二进制文件，获取系数矩阵 fclose(file);%关闭二进制头文件 tic;%计时开始

20

计算方法B上机报告――题目五

%对非压缩带状矩阵求解 if ver=='102'

a=zeros(n,n);

%提取系数矩阵A for i=1:n

for j=1:n

a(i,j)=A((i-1)*n+j); end end

%提取矩阵b b=zeros(n,1); for i=1:n

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

%列主元Gauss消去法 %消去过程 for k=1:n-1 m=k;

%寻找第k步消去时的主元 for i=k+1:n

if abs(a(m,k))

if a(m,k)==0

disp('An error occured while judge the value of a'); return end

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

temp=a(k,j); a(k,j)=a(m,j); a(m,j)=temp; temp=b(k); b(k)=b(m); b(m)=temp; end

for i=k+1:n

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

21

计算方法B上机报告――题目五

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

%对压缩带状矩阵求解

elseifver=='202'%高斯消去法 m=p+q+1; a=zeros(n,m);

%提取系数矩阵A for i=1:1:n for j=1:1:m

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

%提取矩阵b b=zeros(n,1); for i=1:1:n

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

%列主元Gauss消去法 %消去过程 for k=1:1:(n-1)

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

disp('An error occured while judge the value of a'); break; end

length1=n; if (k+p)

length1=k+p; end

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

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 end

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

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

22

计算方法B上机报告――题目五

sum=0;

for k=(n-1):-1:1 sum=b(k); length2=n; if (k+q)

length2=k+q; end

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

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

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

%输出方程组的解

disp('方程组的的解为：') x

toc;%计时结束

5. 计算结果

(1) Dat51.dat 15阶非压缩带状方程组计算结果

程序运行结果如下所示:

文件标识符为id =F1E1D1A0 文件版本号为ver =102 矩阵A的阶数n =15 矩阵A的上带宽q =3 矩阵A的下带宽p =3 方程组的的解为： x =

1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000

1.00001.00001.00001.0000 时间已过 0.000558 秒。

(2) Dat52.dat 20阶压缩带状方程组计算结果

程序运行结果如下所示:

文件标识符为id =F1E1D1A0 文件版本号为ver =202 矩阵A的阶数n =20 矩阵A的上带宽q =5 矩阵A的下带宽p =5 方程组的的解为： x =

1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 1.00001.00001.00001.0000 1.00001.00001.00001.00001.0000

23

计算方法B上机报告――题目五

时间已过 0.001166 秒。

(3) Dat53.dat 2000阶左右非压缩带状方程组计算结果

程序运行结果如下所示:

文件标识符为id =F1E1D1A0 文件版本号为ver =102 矩阵A的阶数n =2160 矩阵A的上带宽q =5 矩阵A的下带宽p =5

由于解中3.1400数目过多无法判断，通过程序：

n=length(x) j=0; for i=1:n

if abs(x(i)-3.14)<10^-5 j=j+1; end end j

对解的数目进行判断得：

num = 2160 j = 2092

即方程组的的解为:2113个3.1400，68个其他解，此处不一一列出。

时间已过 144.156020 秒。

(4) Dat54.dat 4000阶左右非压缩带状方程组计算结果

程序运行结果如下所示:

文件标识符为id =F1E1D1A0 文件版本号为ver =202 矩阵A的阶数n =43240 矩阵A的上带宽q =4 矩阵A的下带宽p =4

由于解中5.2100数目过多无法判断，通过程序：

num=length(x) j=0;

for i=1:1:num

if abs(x(i)-5.21)<10^-5 j=j+1; end end j

对解的数目进行判断得：

num = 43240 j = 43240

即方程组的的解为:43240个5.2100

时间已过 1.841273 秒。

24

计算方法B上机报告――题目五

6. 分析总结

利用列主元Gauss消去法可以很方便的计算高阶线性方程组，对比2000阶左右的非压缩带状方程组和4000阶左右的压缩带状方程组求解所用的时间，可以看出虽然压缩带状方程组的阶数远大于非压缩带状方程组，但其程序运行的总时间却远小于后者。故，在求解高阶方程组时对矩阵做适当处理可大大节省计算时间。

7. 实际问题

等效电路如下图所示，求解各支路电流及输出电压。

1A10Ω I2-30V +7A10Ω I110Ω +I3-UoI420Ω 02A

解：利用节点电压法取参考节点如图所示，其他3个节点电压分别为Un1、Un2、Un3。列节点电压方程：

??11?130?U?U?1???1010?n110n210?????11?130?U?U?1? ??n2?n11010??1010??1?11??U??n1???Un2?2?1020??20对于线性方程组Ax=b整理得:

0??0.2?0.1?A=??0.10.15?0.05??

??0.050.15??0?b???210?2?T采用上述列主元高斯消去法得：x= 40 100 20

得到个节点电压后则各支路电流为：I1=4A I2=-3A I3=4A 输出电压为：Uo=80V

25

本文来源：https://www.bwwdw.com/article/9vc8.html