计算方法B上机报告

计算方法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