线性方程组的数值解法实验报告
- 格式:doc
- 大小:326.00 KB
- 文档页数:9
迭代法解线性方程组数值分析实验报告一、实验目的本次实验旨在深入研究和掌握迭代法求解线性方程组的基本原理和方法,并通过数值实验分析其性能和特点。
具体目标包括:1、理解迭代法的基本思想和迭代公式的推导过程。
2、掌握雅克比(Jacobi)迭代法、高斯赛德尔(GaussSeidel)迭代法和超松弛(SOR)迭代法的算法实现。
3、通过实验比较不同迭代法在求解不同类型线性方程组时的收敛速度和精度。
4、分析迭代法的收敛性条件和影响收敛速度的因素。
二、实验原理1、线性方程组的一般形式对于线性方程组$Ax = b$,其中$A$ 是$n×n$ 的系数矩阵,$x$ 是$n$ 维未知向量,$b$ 是$n$ 维常向量。
2、迭代法的基本思想迭代法是从一个初始向量$x^{(0)}$出发,按照某种迭代公式逐步生成近似解序列$\{x^{(k)}\}$,当迭代次数$k$ 足够大时,$x^{(k)}$逼近方程组的精确解。
3、雅克比迭代法将系数矩阵$A$ 分解为$A = D L U$,其中$D$ 是对角矩阵,$L$ 和$U$ 分别是下三角矩阵和上三角矩阵。
雅克比迭代公式为:$x^{(k+1)}= D^{-1}(b +(L + U)x^{(k)})$。
4、高斯赛德尔迭代法在雅克比迭代法的基础上,每次计算新的分量时立即使用刚得到的最新值,迭代公式为:$x_i^{(k+1)}=(b_i \sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}\sum_{j=i+1}^{n}a_{ij}x_j^{(k)})/a_{ii}$。
5、超松弛迭代法在高斯赛德尔迭代法的基础上引入松弛因子$\omega$,迭代公式为:$x_i^{(k+1)}= x_i^{(k)}+\omega((b_i \sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}\sum_{j=i}^{n}a_{ij}x_j^{(k)})/ a_{ii} x_i^{(k)})$。
线性方程组求解的数值实验一、课题名称:双对角线线性方程组的数值实验二、班级和姓名:09101900 张争(学号 0910190161) 三、实验问题摘要:考虑一种特殊的对角线元素不为零的双对角线线性方程组(以n=7为例)⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡7665544332211d a d a d a d a d a da d ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡7654321x x x x x x x =⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡7654321b b b b b b b 写出解一般的n (奇数)阶方程组的程序(不要用消元过程,因为不用它可以十分方便的解出这个方程组)。
四、过程中的问题及处理方法:解方程组时应当先解出x1和xn,r 然后根据公式求出其他的解,由于矩阵的特殊性,应该从头尾开始解,最后解出21+n x 。
编程序时使用循环语句要分情况讨论。
五、计算公式: 1x =11d b kk k k k d x a b x 11---=(k<n/2)n x =nnd b kk k k k d x a b x 1+-=(k>n/2)21212123212121+--++++--=n n n n n n n d x a x a b x六、结构程序设计输入矩阵阶数→输入各系数→解出前(n-1)/2个解→解出后(n-1)/2个解→解出第(n+1)/2个解七、程序运行结果:这里令所有的系数等于1,实验结果如下:八、主要程序段介绍:解方程函数:由于矩阵具有好的对称性,解的形式也一样,用循环语句。
但是第一个和最后一个解以及中间一个解的形式不同,需要单独解出,程序如下:void jiefangcheng(){int i;x[1]=b[1]/d[1];for(i=2;i<(i+1)/2;i++)x[i]=(b[i]-a[i-1]*x[i-1])/d[i];x[n]=b[n]/d[n];for(i=n-1;i>(i+1)/2;i--)x[i]=(b[i]-a[i]*x[i+1])/d[i];x[(n+1)/2]=(a[(n-1)/2]*x[(n-1)/2]-a[(n+1)/2]*x[(n+3)/2])/d[(n+1)/2];}九、源程序如下://*****************张争******学号0910190161**********#include<iostream.h>#define N 10000int n;int a[N],b[N],d[N],x[N];void input(){int i;cout<<"请输入方程组的阶数:";cin>>n;cout<<"请输入方程组的系数:"<<endl;for(i=1;i<=n;i++){cout<<"d"<<i<<"=";cin>>d[i];cout<<"b"<<i<<"=";cin>>b[i];}for(i=1;i<n;i++){cout<<"a"<<i<<"=";cin>>a[i];}}void jiefangcheng(){int i;x[1]=b[1]/d[1];for(i=2;i<(i+1)/2;i++)x[i]=(b[i]-a[i-1]*x[i-1])/d[i];x[n]=b[n]/d[n];for(i=n-1;i>(i+1)/2;i--)x[i]=(b[i]-a[i]*x[i+1])/d[i];x[(n+1)/2]=(a[(n-1)/2]*x[(n-1)/2]-a[(n+1)/2]*x[(n+3)/2])/d[(n+1)/2];}void output(){ int i;cout<<"方程组的解为:"<<endl;for(i=1;i<=n;i++)cout<<"x"<<i<<"="<<x[i]<<endl;}void main(){input();jiefangcheng();output();}。
实验报告——线性方程组的数值解法姓名:班级:学号:日期:一 实践目的1. 熟悉求解线性方程组的有关理论和方法。
2. 会编列主元消去法,全主元消去法,雅克比迭代法及高斯-赛德尔迭代法的程序。
3.通过实际计算,进一步了解各种方法的优缺点,选择合适的数值方法。
4. 进一步应用数学知识,扩展数学思维,提高编程能力。
二 问题定义及题目分析1.求解线性方程是实际中常遇到的问题。
而一般求解方法有两种,一个是直接求解法,一个是迭代法。
2.直接法是就是通过有限步四则运算求的方程准确解的方法。
但实际计算中必然存在舍入误差,因此这种方法只能得到近似解。
3.迭代法是先给一个解的初始近似值,然后按一定的法则求出更准确的解,即是用某种极限过程逐步逼近准确解的方法。
4.这次实践,将采用直接法:列主元高斯消去法,全主元高斯消去法;迭代法:雅克比迭代法和高斯-赛德尔迭代法。
三 详细设计 设有线性方程组11112211n n a x a x a x b +++=21122222n n a x a x a x b +++= 1122n n nn n n a x a x a x b +++=令A= 111212122212n n n n nn a a a aa a a a a ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ x=12n x x x ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ b= 12n b b b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦一. 上三角方程组的求解很简单。
所以若能将方程组化为上三角方程组,那就很容易得到方程的解,高斯消去法则是一种简洁高效的方法,可以将方程组化为上三角方程组,但是这种方法必须满足每一步时0kk a =才能求解,还有当kk a 绝对值很小时,作分母会引起较大的舍入误差。
因此在消元过程中应尽量选择绝对值较大的系数作为主元素。
1.列主元消去法很好的解决这一问题。
将1x 的系数1(1)k a k n ≤≤中绝对值最大者作为主元素,交换第一行和此元素所在的行,然后主元消去非主元项1x 的系数,。
1大学数学实验 实验报告 | 2014/4/5一、 实验目的1、学习用Matlab 软件数值求解线性代数方程组,对迭代法的收敛性和解的稳定性作初步分析;2、通过实例学习用线性代数方程组解决简化问题。
二、 实验内容项目一:种群的繁殖与稳定收获:种群的数量因繁殖而增加,因自然死亡而减少,对于人工饲养的种群(比如家畜)而言,为了保证稳定的收获,各个年龄的种群数量应维持不变。
种群因雌性个体的繁殖而改变,为方便起见以下种群数量均指其中的雌性。
种群年龄记作k=1,2,…,n ,当年年龄k 的种群数量记作x k ,繁殖率记作b k (每个雌性个体1年的繁殖的数量),自然存活率记作s k (s k =1−d k ,d k 为1年的死亡率),收获量记作ℎk ,则来年年龄k 的种群数量x ̌k 应该为x ̌k =∑b k n k=1x k , x ̌k+1=s k x k −ℎk , (k=1,2,…,n -1)。
要求各个年龄的种群数量每年维持不变就是要求使得x ̌k =x k , (k=1,2,…,n -1).(1) 如果b k , s k 已知,给定收获量ℎk ,建立求各个年龄的稳定种群数量x k 的模型(用矩阵、向量表示).(2) 设n =5,b 1=b 2=b 5=0,b 3=5,b 4=3,s 1=s 4=0.4,s 2=s 3=0.6,如要求ℎ1~ℎ5为500,400,200,100,100,求x 1~x 5.(3) 要使ℎ1~ℎ5均为500,如何达到?问题分析:该问题属于简单的种群数量增长模型,在一定的条件(存活率,繁殖率等)下为使各年龄阶段的种群数量保持不变,各个年龄段的种群数量将会满足一定的要求,只要找到种群数量与各个参量之间的关系,建立起种群数量恒定的方程就可以求解出各年龄阶段的种群数量。
模型建立:根据题目中的信息,令x ̌k =x k ,得到方程组如下:{x ̌1=∑b k nk=1x k =x 1x ̌k+1=s k x k −ℎk =x k+1整理得到:{−x 1∑b k nk=1x k =0−x k+1+s k x k =ℎk2 大学数学实验 实验报告 | 2014/4/52写成系数矩阵的形式如下:A =[b 1−1b 2b 3s 1−100s 2−1…b n−1b n0000⋮⋱⋮000000000⋯00−10s n−1−1]令h =[0, ℎ1,ℎ2,ℎ3,…,ℎn−2,ℎn−1]Tx =[x n , x n−1,…,x 1]T则方程组化为矩阵形式:Ax =h ,即为所求模型。
第12 章 实验十一线性方程组的数值解法实验目的:理解线性方程组基本定理,掌握求解线性方程组的各种方法,逆矩阵求法、高斯消去法、带选主元的分解法、迭代法。
12.1 线性方程组(基本定理) 已知m 个方程n 个未知量的方程组:⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++m n mn m m n n n n b x a x a x a b x a x a x a b x a x a x a 22112222212********* (12.1)其中ij a 为已知系数,j x 未知,j b 为已知项,当),,2,1(n j b j =全为零时,称(12.1)方程组为线性齐次方程组,当),,2,1(n j b j =不全为零时,称(12.1)方程组为线性非齐次方程组。
上述方程组也可表示如下矩阵形式:b Ax = 其中A 、x 和y 分别定义为:⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=mn m n a a a a A LL L L L 1111⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=n x x x 1 ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=m b b b 1。
方程组(12.1)中的m 与n 有三种情况:m=n 或m<n 或m>n ,m=n 是最常见的。
当m=n 时,有系数矩阵的行列式0≠A ,方程组(12.1)有唯一解。
用矩阵的秩可以得到更满意的结果:(1)当n B R A R ==)()(时,方程组(12.1)有唯一解, (2)当n B R A R <=)()(时,方程组(12.1)有无穷多组解, (3)当)()(B R A R ≠时,方程组(12.1)无解。
12.2举例例12.1 用MATLAB 求解线性方程组:⎪⎪⎩⎪⎪⎨⎧=+++-=+++=+++=+++62332022428340213424321432143214321x x x x x x x x x x x x x x x x解法一:逆矩阵求法设 A=[1 2 1 4 ; 2 0 4 3 ;4 2 2 1 ;-3 1 3 2 ]; det(A)=-180, B=[13 28 20 6]’; 则:X=A\B % A\表示矩阵A 的逆,或用1-A 表示得到X =3.0000 -1.00004.0000 2.0000也可以用如下命令: X=B ’/A ’ 得到X =3.0000 -1.00004.0000 2.0000解法二:在MA TLAB 上,用高斯消去法求解方程组:⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-62028132313122434024214321x x x x (12.2) 首先写出方程组对应的增广矩阵:a=[ 1 2 1 4 13; ... 2 0 4 3 28 ; ... 4 2 2 1 20 ; ...-3 1 3 2 6 ]其中前四列是系数矩阵,最后一列是方程组的常数列。
线性代数方程组的数值解法【实验目的】1. 学会用MATLAB 软件数值求解线性代数方程组,对迭代法的收敛性和解的稳定性作初步分析;2. 通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】【题目1】通过求解线性方程组A1x=b1和A2x=b2,理解条件数的意义和方程组的性态对解的影响。
其中A1是n阶范德蒙矩阵,即⎡1x0⎢1x1⎢A1=⎢⎢⎢⎣1xn-12x0x12 2xn-1n-1⎤ x0⎥ x1n-1⎥1,...,n-1 ,xk=1+0.1k,k=0,⎥ n-1⎥ xn-1⎥⎦A2是n阶希尔伯特矩阵,b1,b2分别是A1,A2的行和。
(1)编程构造A1(A2可直接用命令产生)和b1,b2;你能预先知道方程组A1x=和A2x=。
b2的解吗?令n=5,用左除命令求解(用预先知道的解可检验程序)b1(2)令n=5,7,9,…,计算A1,A2的条件数。
为观察它们是否病态,做以下试验:b1,b2不变,A1和A2的元素A1(n,n),A2(n,n)分别加扰动ε后求解;A1和A2不变,b1,b2的分量b1(n),分析A和b的微小扰动对解的影响。
b2(n)分别加扰动ε求解。
ε取10-1010,-8,10-6。
(3)经扰动得到的解记做x~,计算误差-x~x,与用条件数估计的误差相比较。
1.1构造A1,A2和b1,b2首先令n=5,构造出A1,A2和b1,b2。
首先运行以下程序,输出A1。
运行以下程序对A1,A2求行和:由于b1,b2分别是A1,A2的行和,所以可以预知x1=运行下列程序,用左除命令对b1,b2进行求解:得到以下结果: T。
x2=(1,1, ,1)1.2 计算条件数并观察是否为病态1.不加扰动,计算条件数。
运行以下程序:由此可知,A1,A2的条件数分别是3.574∗10, 4,766∗10。
2.b1,b2不变,A1(n,n),A2(n,n)分别加扰动(1)n=5时设x11,x12,x13分别为A1添加扰动10−10,10−8,10−6后的解。
线性方程组的数值解法一 实验目的1 学会用MATLAB 软件求解线性代数方程组,对迭代法的收敛性和解的稳定性做初步分析;2 通过实例学习线性代数方程组解决简化的实际问题。
二 实验内容1 已知方程组Ax =b ,其中A ∈ℝ20×20,定义为:31/21/41/231/21/41/41/231/21/41/41/231/21/41/23--⎡⎤⎢⎥---⎢⎥⎢⎥---⎢⎥-⎢⎥⎢⎥---⎢⎥--⎣⎦试通过迭代法求解此方程组,认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对收敛速度的影响。
实验要求:(1) 选取不同的初始向量x (0)和不同的方程组右端项向量b ,给定迭代误差要求,用雅可比迭代法和高斯-赛德尔迭代法计算,观测得到的迭代向量序列是否均收敛?若收敛,记录迭代次数,分析计算结果并得出你的结论。
(2) 取定右端项向量b 和初始向量x (0),将A 的主对角线元素成倍增长若干次,非对角线元素不变,每次用雅可比迭代法计算,要求迭代误差满足‖x (k+1)−x (k)‖∞<10−5,比较收敛速度,分析现象并得出你的结论。
初步解决:首先建立利用雅克比迭代法和高斯-赛德尔迭代法计算的函数M 文件。
雅克比迭代法:高斯-赛德尔迭代法:关于迭代法的收敛性:由原式Ax=b化简可得x=Bx+f,也就是说,研究此处迭代法是否收敛实际上就是研究有原矩阵A经变换之后所得矩阵B的性质。
而由于该题第一问中,A保持不变,所以对于收敛性只需研究一次。
在命令栏中输入以下命令:(1)构造题目中的矩阵A:(2)提取对角矩阵、上三角矩阵、下三角矩阵:(3)对于两种迭代方法分别求B:(4)分别计算两个矩阵的三种范数:n1和m1、n2和m2、n3和m3分别是矩阵B1、B2的1-范数、2-范数和∞范数。
由书上定理,矩阵的谱半径不超过任何一个范数,即ρ(B)≤‖B‖,而由图中可以看出,两个矩阵的六个范数没有一个大于1,所以两个矩阵的谱半径一定都小于1,所以此时两种迭代方法均收敛。
数值分析实验报告(一)一、实验名称:求解线性方程组二、实验问题:从2到9的每一个n 值,解n 阶方程组Ax=b 。
这里A 和b 如下定义:a ij =(i+j-1)-1,b i =p(n+i-1)-p(i-1)式中p(x)=(x 2/24)(2+x 2(-7+n 2(14+n(12+3n))))解释发生的现象。
三、计算公式(数学模型):Gauss 消去法是将方程组Ax=b ,用逐次消去未知数的方法,化为等价的三角方程组,再回代求解的方法。
设Ax=b ,A n n R ⨯∈,如果),...,2,1(,0)(n k a k kk=≠,则可通过Gauss 消去法将Ax=b 约化为等价的三角形方程组⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡)()2(2)1(121)()2(2)2(22)1(1)1(12)1(11n n n n nn n n b b b x x x a a a a a a且计算公式为:(a)消元计算(k=1,2,…,n)()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧+=-=+=-=+==++n k i a m b b n k j i a m a a n k i a a m k k ik k i k i k kj ik k ij k ijk kn k ik ik ,,1,,1,,1)()()1()()()1()()( (b)回代计算:()⎪⎪⎩⎪⎪⎨⎧-=-==∑+=n i j i ii ji i ij i i in nn n n n n ia x ab x a b x 1)()()()()(1,2,,1四、结构程序设计:1、数据流向图2五、计算程序C++程序//改进的平方根法#include<iostream.h>double p(double x,double n){return x*x/24*(2+x*x*(-7+n*n*(14+n*(12+3*n))));}void main(void){cout<<"*************************线性方程组求解的数值实验题*****************************"<<'\n';double A[9][9]={0};double L[9][9]={0};double B[9]={0};double d[9]={0};double x[9]={0};double y[9]={0};double dd=0;int n,i,j,k;for(n=2;n<10;n++){cout<<n<<"阶方程组Ax=b的解为"<<'\n';for(i=0;i<n;i++){for(j=0;j<n;j++){A[i][j]=1.0/(i+j+1); //赋Ax=b中系数A的初值(A是对称阵)}B[i]=p((double)n+i,(double)n)-p((double)i,(double)n); //赋Ax=b中系数b的初值(调用函数p)}d[0]=A[0][0]; //d[]用来记录A[k][k]for(i=0;i<n;i++) L[i][0]=A[i][0]/d[0]; //LDL'分解中L[i][j]的第一列for(i=1;i<n;i++){for(k=0;k<i;k++) dd+=L[i][k]*L[i][k]*d[k];d[i]=A[i][i]-dd; //LDL'分解中的d[i]dd=0;for(j=i;j<n;j++){for(k=0;k<j;k++) dd+=L[i][k]*d[k]*L[j][k]; //LDL'分解中除第一列外的L[i][j]L[i][j]=(A[i][j]-dd)/d[j];dd=0;}}y[0]=B[0]; //对y的回代过程for(i=1;i<n;i++){for(k=0;k<i;k++) dd+=L[i][k]*y[k];y[i]=B[i]-dd;dd=0;}x[n-1]=y[n-1]/d[n-1]; //对x的回代过程for(i=n-2;i>=0;i--){for(k=i+1;k<n;k++) dd+=L[k][i]*x[k];x[i]=y[i]/d[i]-dd;dd=0;}for(i=0;i<n;i++){cout<<'x'<<i+1<<"="<<x[i]<<'\n';}cout<<'\n';}} //在C++数组中是以开始的,所以i,j,k取值相应变到从0开始运行结果如下:六、结果和讨论:对于2到9阶的线性方程组的解来看,Gauss 消去法有其自己的弊端,使有些结果有误差。
实验5 线性代数方程组的数值解法化工系 分0班 2010011805 张亚清【实验目的】1、 学会用MATLAB 软件数值求解线性代数方程组,对迭代法的收敛性和解的稳定性作初步分析;2、 通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】 1、 题目3已知方程组Ax=b ,其中2020⨯ℜ∈A ,定义为试通过迭代法求解此方程组,认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对收敛速度的影响。
实验要求:(1)选取不同的初始向量)0(x和不同的方程组右端向量b ,给出迭代误差要求,用雅克比迭代法和高斯-赛德尔迭代法计算,观测得到的迭代向量序列是否均收敛?若收敛,记录迭代次数,分析计算结果并得出你的结论;(2)取定右端向量b 和初始向量)0(x,将A 的主对角线元素成倍增长若干次,非主对角线元素不变,每次用雅克比迭代法计算,要求迭代误差满足5)()1(10-∞+<-k k x x ,比较收敛速度,分析现象并得出你的结论。
【问题分析】对于线性方程组Ax=b ,满足0≠ii a (i=1,2,...,n ),将A 分解为A=D-L-U ,则雅克比迭代公式等价于如下的矩阵形式,...2,1,0,)(1)(1)1(=++=--+k b D x D L D x k k或b D f U L D B f x B xJ J J k J k 11)()1(),(,--+=+=+=。
类似的,Ax=b 的高斯-赛德尔迭代公式等价于如下矩阵形式b L D f U L D B f x B x J S G S G k S G k 11)()1()(,)(,-----+-=-=+=。
【问题解答】(1)选取初始向量T x )1,...,1,1()0(=,T b )1,...,1,1(=,迭代要求为4)()1(10-∞+<-k k x x 。
将A 按A=D-L-U 分解为如下三个矩阵:①对方程组进行雅克比迭代,利用MATLAB 编程计算。
实验三线性代数方程组的数值解法一、迭代法求解方程组㈠问题描述给定方程组的矩阵A,通过迭代法求解方程组。
1、选取不同的初始向量和不同的右端项向量,给定误差要求,用两种迭代法计算;2、去顶右端项向量和初始向量,将A的主对角线元素成倍增长若干次,非主对角线元素不变,用雅克比迭代法计算。
㈡方法与公式1、雅克比迭代法2、高斯-赛德尔迭代法㈢结果与分析1、不同初始向量、不同右端项向量、不同精度要求(1)初始向量定为zeros(n,1);①b=zeros(n,1)迭代次数为0,直接得到结果。
③b = 1:n(2)初始向量定为one s(n,1)①b=zeros(n,1)事实上,迭代100次时,所得结果约为10^-32,已经可以认为是0,但是由于没有达到精度要求,故不算收敛。
②b=ones(n,1)④b = n:1(3)初始向量定为1:n①b=zeros(n,1)②b=ones(n,1)(4)初始向量定为n:1①b=zeros(n,1)②b=ones(n,1)④b = n:1(5)简要小结a.在个别情况下雅可比迭代法收敛速度极慢,但事实上没有达到收敛时其计算结果已经可以接受;b.要求的精度越高,迭代的次数越多,迭代的次数与所要求的精度的对数值近似呈线性,也就是说两者近似呈指数关系;b.高斯-赛德尔迭代法有着比雅可比更好的迭代特性;2、更改A的主对角线元素(1) b =20:1;初值 1:20(2) b =20:1;初值 20:1(3) b =1:20;初值 20:1(4) b =[3;5;2;6;8;23;5;8;32;63;23;5;2;12;0;23;1;564;2;65]; 初值 ones(20,1)(5)简要小结a.迭代的次数随着对角线元素的成倍的增长而降低,趋于一稳定值;b.右端项以及迭代初值仅当对角线元素较小时对迭代次数起有作用,对角线元素成倍数增加后,迭代次数不变。
3、总结由以上各个对比可以得出以下结论:a.使用迭代法求解方程组时时,要求的精度越高,迭代次数越大;b.高斯-赛德尔迭代法的迭代次数要比雅可比迭代法迭代次数低;c.雅可比迭代的次数随着矩阵A对角线元素的成倍的增长而降低,d.当矩阵A的对角线元素足够大时,雅可比迭代法的迭代次数趋于稳定值;㈣程序清单1、第一问中的雅可比迭代function [y,k] = jacobi(A,b,m,tol)D = diag(diag(A));L = - tril(A,-1);U = - triu(A,1);n = length(A);y = ones(n,1);BJ=D\(L+U);fJ=D\b;k=0;while norm(A*y-b)/norm(b)>tol && k<mk = k+1;y = BJ*y+fJ;end2、高斯-赛德尔迭代法function [y,k] = GuassSeidel(A,b,m,tol)D = diag(diag(A));L = - tril(A,-1);U = - triu(A,1);n = length(A);BG = (D-L)\U;FG = (D-L)\b;y = zeros(n,1);k=0;while norm(A*y-b)/norm(b)>tol && k<mk = k+1;y = BG*y+FG;end3、第二问中的雅可比迭代function [y,k] = jacobi1(A,b,m,tol) D = diag(diag(A));L = - tril(A,-1);U = - triu(A,1);n = length(A);yk = ones(20,1);BJ=D\(L+U);fJ=D\b;k=0;yk1 = BJ*yk+fJ;k = k+1;ttol = norm(yk1-yk,inf);while ttol>tol && k<mk = k+1;yk = yk1;yk1 = BJ*yk+fJ;ttol = norm(yk1-yk,inf);endy=yk1;4、第一问脚本n=20;A1 = sparse(1:n,1:n,3,n,n);A2 = sparse(1:n-1,2:n,-1/2,n,n);A3 = sparse(2:n,1:n-1,-1/2,n,n);A4 = sparse(1:n-2,3:n,-1/4,n,n);A5 = sparse(3:n,1:n-2,-1/4,n,n);A = A1+A2+A3+A4+A5;b =zeros(20,1);m=10000;for i = 3:10tol = 10^(-i);[xJ,k1(i)] = jacobi(A,b,m,tol); [xG,k2(i)] = GuassSeidel(A,b,m,tol); endk1k25、第二问脚本n=20;k = 1;A1 = sparse(1:n,1:n,3,n,n);A2 = sparse(1:n-1,2:n,-1/2,n,n);A3 = sparse(2:n,1:n-1,-1/2,n,n);A4 = sparse(1:n-2,3:n,-1/4,n,n);A5 = sparse(3:n,1:n-2,-1/4,n,n);b =[3;5;2;6;8;23;5;8;32;63;23;5;2;12;0;23;1;564;2;65]; m=10000;tol = 10^(-5);for k = 1:10A = k*A1+A2+A3+A4+A5;[xJ,k1(k)] = jacobi1(A,b,m,tol);endk1二、一年生植物的繁殖㈠问题描述对于5.1.2的假设,给定参数,试分别用追赶法、稀疏系数矩阵和满矩阵求解;若b 有10%误差,估计对结果的影响。
重庆大学学生实验报告实验课程名称计算方法开课实验室DS1421学院年级专业学生姓名学号开课时间至学年第学期1.实验目的(1)高斯列主元消去法求解线性方程组的过程(2)熟悉用迭代法求解线性方程组的过程(3)设计出相应的算法,编制相应的函数子程序2.实验内容分别用高斯列主元消去法 ,Jacobi 迭代法,Gauss--Saidel 迭代法,超松弛迭代法求解线性方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------725101391444321131243301024321x x x x 3.实验过程解:(1)高斯列主元消去法编制高斯列主元消去法的M 文件程序如下:%高斯列主元消元法求解线性方程组Ax=b%A 为输入矩阵系数,b 为方程组右端系数%方程组的解保存在x 变量中format long;%设置为长格式显示,显示15位小数A=[2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13]b=[10,5,-2,7]'[m,n]=size(A);%先检查系数正确性if m~=nerror('矩阵A 的行数和列数必须相同');return;endif m~=size(b)error('b 的大小必须和A 的行数或A 的列数相同');return;end%再检查方程是否存在唯一解if rank(A)~=rank([A,b])error('A 矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');return;endc=n+1;A(:,c)=b; %(增广)for k=1:n-1[r,m]=max(abs(A(k:n,k))); %选主元m=m+k-1; %修正操作行的值if(A(m,k)~=0)if(m~=k)A([k m],:)=A([m k],:); %换行endA(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); %消去 endendx=zeros(length(b),1); %回代求解x(n)=A(n,c)/A(n,n);for k=n-1:-1:1x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k);enddisp('X=');disp(x);format short;%设置为默认格式显示,显示5位运行,结果如下所示:(2)Jacobi迭代法编制迭代计算的M文件程序如下:%Jacobi迭代法求解% A为方程组的增广矩阵clc;A=[2,10,0,-3,10;-3,-4,-12,13,5;1,2,3,-4,-2;4,14,9,-13,7] MAXTIME=50;%最多进行50次迭代eps=1e-5;%迭代误差[n,m]=size(A);x=zeros(n,1);%迭代初值y=zeros(n,1);k=0;%进入迭代计算disp('迭代过程X的值情况如下:')disp('X=');while 1disp(x');for i=1:1:ns=0.0;for j=1:1:nif j~=is=s+A(i,j)*x(j);endy(i)=(A(i,n+1)-s)/A(i,i);endendfor i=1:1:nmaxeps=max(0,abs(x(i)-y(i))); %检查是否满足迭代精度要求 endif maxeps<=eps%小于迭代精度退出迭代for i=1:1:nx(i)=y(i);%将结果赋给xendreturn;endfor i=1:1:n%若不满足迭代精度要求继续进行迭代x(i)=y(i);y(i)=0.0;endk=k+1;if k>MAXTIME%超过最大迭代次数退出error('超过最大迭代次数,退出');return;endend运行该程序结果如下:(3)Gauss--Saidel迭代法编制求解程序Gauss_Seidel.m如下:%Gauss_Seidel.m% A为方程组的增广矩阵clc;format long;A=[2,10,0,-3,10;-3,-4,-12,13,5;1,2,3,-4,-2;4,14,9,-13,7][n,m]=size(A);%最多进行50次迭代Maxtime=50;%控制误差Eps=10E-5;%初始迭代值x=zeros(1,n);disp('x=');%迭代次数小于最大迭代次数,进入迭代for k=1:Maxtimedisp(x);for i=1:ns=0.0;for j=1:nif i~=js=s+A(i,j)*x(j);%计算和endendx(i)=(A(i,n+1)-s)/A(i,i);%求出此时迭代的值end%因为方程的精确解为整数,所以这里将迭代结果向整数靠近的误差作为判断迭代是否停止的条件if sum((x-floor(x)).^2)<Epsbreak;end;end;X=x;disp('迭代结果:');Xformat short;运行结果如下所示:(4)超松弛迭代法编写函数M文件如下:%SOR法求解%w=1.45%方程组系数矩阵clc;A=[2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13]%方程组右端系数b=[10,5,-2,7]'w=1.45;%最大迭代次数Maxtime=100;%精度要求Eps=1E-5;%以15位小数显示format long;n=length(A);k=0;%初始迭代值x=ones(n,1);y=x;disp('迭代过程:');disp('x=');while 1y=x;disp(x');%计算过程for i=1:ns=b(i);for j=1:nif j~=is=s-A(i,j)*x(j);endendif abs(A(i,i))<1E-10 | k>=Maxtimeerror('已达最大迭代次数或矩阵系数近似为0,无法进行迭代'); return;ends=s/A(i,i);x(i)=(1-w)*x(i)+w*s;endif norm(y-x,inf)<Eps%达到精度要求退出计算break;endk=k+1;enddisp('最后迭代结果:');%最后的结果X=x'%设为默认显示格式format short; 结果如下:4.实验环境及实验文件存档名实验环境:Matlab7.0文件存档名:Gauss.m,Jacobi.m,Gauss_Seidel.m,SOR.m5.实验结果及分析=1.0000=2.0000=3.0000=4.0000经过验证,高斯列主元消结果正确。
实验5 线性代数方程组的数值解法分1 黄浩 43一、实验目的1.学会用MATLAB软件数值求解线性代数方程组,对迭代法的收敛性和解的稳定性作初步分析;2.通过实例学习用线性代数方程组解决简化的实际问题。
二、For personal use only in study and research; not forcommercial use三、四、实验内容1.《数学实验》第二版(问题1)问题叙述:通过求解线性方程组,理解条件数的意义和方程组性态对解的影响,其中是n阶范德蒙矩阵,即是n阶希尔伯特矩阵,b1,b2分别是的行和。
(1)编程构造(可直接用命令产生)和b1,b2;你能预先知道方程组和的解吗?令n=5,用左除命令求解(用预先知道的解可验证程序)。
(2)令n=5,7,9,…,计算和的条件数。
为观察他们是否病态,做以下试验:b1,b2不变,和的元素,分别加扰动后求解;和不变,b1,b2的分量b1(n),b2(n)分别加扰动后求解。
分析A与b的微小扰动对解的影响。
取10^-10,10^-8,10^-6。
(3)经扰动得到的解记做,计算误差,与用条件数估计的误差相比较。
模型转换及实验过程:(1)小题.由b1,b2为,的行和,可知方程组和的精确解均为n 行全1的列向量。
在n=5的情况下,用matlab编程(程序见四.1),构造,和b1,b2,使用高斯消去法得到的解x1,x2及其相对误差e1,e2(使用excel计算而得)为:由上表可见,当n=5时,所得的解都接近真值,误差在10^-12的量级左右。
(2)小题分别取n=5,7,9,11,13,15,计算和的条件数c1和c2,(程序见四.2),结果如下:由上表可见,二者的条件数都比较大,可能是病态的。
为证实和是否为病态,先保持b不变,对做扰动,得到该情况下的高斯消元解,(程序见四.3),结果如下:(为使结果清晰简洁,在此仅列出n=5,9,13的情况,n=7,11,15略去)=10^-10时:=10^-8时:=10^-6时:由上表可见:a)对于希尔伯特阵,随着阶数的增加,微小扰动对解带来的影响越来越大,到了n=9时,已经有了6倍误差的解,到了n=13时,甚至出现了22倍误差的解元素;而随着的增加,解的偏差似乎也有增加的趋势,但仅凭上述表格无法具体判断(在下一小题中具体叙述)。
实验5:线性方程组的数值解法化学工程系分2 安振华2012011837【实验目的】1、掌握线性方程组的常用数值解法,包括高斯消去法、LU分解法以及校正法。
2、体验数值计算的时间复杂度和计算规模的关系。
3、加深对数值计算误差的理解。
4、学习使用迭代法等算法,求解非线性方程。
5、学习如何使用MATLAB解非线性方程组和方程组。
【实验容】【实验五:习题9】种群的繁殖与稳定收获:种群的数量因繁殖而增加,因自然死亡而减少,对于人工饲养的种群(比如家畜)而言,为了保证稳定的收获,各个年龄的种群数量应保持不变,种群因雌性个体的繁殖而改变,为方便起见以下种群数量均指其中的雌性。
种群年龄记作k=1,2,…,n,当年年龄k的种群数量记作x k,繁殖率记作b k(每个雌性个体在1年繁殖的数量),自然存活率记作s k(s k=1-d k,d k为1年的死亡率),收获量记作h k,则来年年龄k的种群数量k x应为:111,(1,2,,1)n k k k k k k k x b x x s x h k n +===-=⋅⋅⋅-∑要求各个年龄的种群数量每年维持不变就是要使(1,2,,)k k x x k n ==⋅⋅⋅(1) 若b k ,s k 已知,给定收获量h k ,建立求各年龄的稳定种群数量x k 的模型(用矩阵向量表示)(2) 设n=5,b 1=b 2=b 5=0,b 3=5,b 4=3,s 1=s 4=0.4,s 2=s 3=0.6,如果要求h 1~h 5为500,400,200,100,100,求x 1~x 5 (3) 要使h 1~h 5均为500,如何达到? 【分析】为方便起见以下种群数量均指其中的雌性。
我们并且有以下的假设:(1)雌性个体的繁殖率和存活率在特定的时间是不变的。
(2)人工饲养的种群在质量和数量上是不受外界环境和资源的限制的。
(3)模型中不考虑人为的或是自然的灾害所造成的种群数量、繁殖率和存活率的变动。
姓名:学号:班级:线性方程组的数值解法实习目的(1)通过实习进一步掌握高斯消去法、列主元高斯消去法、柯朗分解法、追赶法以及雅克比迭代法和高斯-赛德尔迭代法的基本思想;(2)通过实习进一步掌握高斯消去法、列主元高斯消去法、柯朗分解法、追赶法以及雅克比迭代法和高斯-赛德尔迭代法的计算步骤,并能灵活应用;(3)通过对高斯消去法、列主元高斯消去法、柯朗分解法、追赶法以及雅克比迭代法和高斯-赛德尔迭代法的调试练习,进一步体会各种算法的特点;(4)通过对上机调试运行,组不培养解决实际问题的编程能力。
实习要求(1)熟悉Turbo C的编译环境;(2)实习前复习高斯消去法、列主元高斯消去法、直接三角分解法、雅克比迭代法、高斯-赛德尔迭代法以及追赶法的基本思想和过程;(3)实习前复习高斯消去法、列主元高斯消去法、直接三角分解法、雅克比迭代法、高斯-赛德尔迭代法以及追赶法的计算步骤。
实习设备(1)硬件设备:单机或网络环境下的微型计算机一台;(2)软件设备:DOS 3.3以上操作系统,Turbo C 2.0编译器。
实习内容实习一高斯消去法(1)用高斯消去法求解线性方程组:(2)要求:①将线性方程组写成用矩阵表示的形式,即Ax=b的形式。
②输出系数矩阵的原始元素和经高斯消去法消去后的矩阵元素。
③经高斯消去法消去后的矩阵是一个什么形式的矩阵?④请写出程序的运行结果。
#include <stdio.h>#define N 4int Gauss(float a[N][N],float b[N]){int i,j,k,flag=1;float t;for(i=0;i<N-1;i++){if(a[i][i]==0){flag=0;break;}else{for(j=i+1;j<N;j++){t=-a[j][i]/a[i][i];b[j]=b[j]+t*b[i];for(k=i;k<N;k++)a[j][k]=a[j][k]+t*a[i][k];}}}return(flag);}void zg_matric(float a[N][N],float b[N]){int i,j;for(i=0;i<N;i++){for(j=0;j<N;j++)printf("%10f",a[i][j]);printf("%10f",b[i]);printf("\n");}printf("\n");}void main(){static float a[N][N]={{1,-4,6,-5},{-5,21,-33,32},{6,-26,43,-48},{5,-24,45,-64}};float b[N]={8,-32,25,-10};float x[N]={0,0,0,0};int i,j,flag;zg_matric(a,b);flag=Gauss(a,b);if(flag==0)printf("Gauss method does not run.");else{x[N-1]=b[N-1]/a[N-1][N-1];for(i=N-2;i>=0;i--){x[i]=b[i];for(j=i+1;j<N;j++)x[i]=x[i]-a[i][j]*x[i];x[i]=x[i]/a[i][i];}for(i=0;i<N;i++)printf("x%d=%11.7f\n",i+1,x[i]);}}(3)思考题:①高斯消去法的基本思想是什么?答:高斯消去法是最古老的求解线性代数方程组的方法之一,是消去法的一种特殊形式。
线性代数方程组的数值解法实验1. 主元的选取与算法的稳定性问题提出:Gauss 消去法是我们在线性代数中已经熟悉的。
但由于计算机的数值运算是在一个有限的浮点数集合上进行的,如何才能确保Gauss 消去法作为数值算法的稳定性呢?Gauss 消去法从理论算法到数值算法,其关键是主元的选择。
主元的选择从数学理论上看起来平凡,它却是数值分析中十分典型的问题。
实验内容:考虑线性方程组 n n n R b R A b Ax ∈∈=⨯,,编制一个能自动选取主元,又能手动选取主元的求解线性方程组的Gauss 消去过程。
实验要求:(1)取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ,则方程有解T x )1,,1,1(* =。
取n=10计算矩阵的条件数。
让程序自动选取主元,结果如何?(2)现选择程序中手动选取主元的功能。
每步消去过程总选取按模最小或按模尽可能小的元素作为主元,观察并记录计算结果。
若每步消去过程总选取按模最大的元素作为主元,结果又如何?分析实验的结果。
(3)取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异,说明主元素的选取在消去过程中的作用。
(4)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。
重复上述实验,观察记录并分析实验结果。
1.1程序清单n=input('矩阵A 的阶数:n=');A=6*diag(ones(1,n))+diag(ones(1,n-1),1)+8*diag(ones(1,n-1),-1); b=A*ones(n,1);p=input('计算条件数使用p-范数,p='); cond_A=cond(A,p) [m,n]=size(A); Ab=[A b];r=input('选主元方式(0:自动;1:手动),r=');Abfor i=1:n-1switch rcase(0)[aii,ip]=max(abs(Ab(i:n,i)));ip=ip+i-1;case (1)ip=input(['第',num2str(i),'步消元,请输入第',num2str(i),'列所选元素所处的行数:']);end;Ab([i ip],:)=Ab([ip i],:);aii=Ab(i,i);for k=i+1:nAb(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i:n+1);end;if r==1Abendend;x=zeros(n,1);x(n)=Ab(n,n+1)/Ab(n,n);for i=n-1:-1:1x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);endx1.2运行结果(1)n=10,矩阵的条件数及自动选主元Cond(A,1) =2.5575×103Cond(A,2) = 1.7276×103Cond(A,inf) =2.5575×103程序自动选择主元(列主元)a.输入数据矩阵A的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=0b.计算结果x=[1,1,1,1,1,1,1,1,1,1]T (2)n=10,手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=10 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:2 …(实际选择时,第k 步选择主元处于第k 行) 最终计算得x=[1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000001, 0.999999999999998, 1.000000000000004, 0.999999999999993, 1.000000000000012, 0.999999999999979, 1.000000000000028]Tb. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=10 计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:3 …(实际选择时,第k 步选择主元处于第k+1行) 最终计算得x=[1,1,1,1,1,1,1,1,1,1]T (3)n=20,手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=20 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:2 …(实际选择时,第k 步选择主元处于第k 行) 最终计算得x=[1.000000000000000,1.000000000000000,1.000000000000000,1.000000000000001,0.999999999999998,1.000000000000004,0.999999999999993,1.000000000000014,0.999999999999972,1.000000000000057,0.999999999999886,1.000000000000227,0.999999999999547,1.000000000000902,0.999999999998209,1.000000000003524,0.999999999993179,1.000000000012732,0.999999999978173,1.000000000029102]T b. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=20 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:3…(实际选择时,第k步选择主元处于第k+1行)最终计算得x=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]T(4)A分别为幻方矩阵,Hilbert矩阵,pascal矩阵和随机矩阵将计算结果列于下表:1.3简要分析计算(1)表明:对于同一矩阵,不同范数定义的条件数是不同的;Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。
实验报告
——线性方程组的数值解法
姓名:
班级:
学号:
日期:
一 实践目的
1. 熟悉求解线性方程组的有关理论和方法。
2. 会编列主元消去法,全主元消去法,雅克比迭代法及高斯-赛德尔迭代法的程序。
3.通过实际计算,进一步了解各种方法的优缺点,选择合适的数值方法。
4. 进一步应用数学知识,扩展数学思维,提高编程能力。
二 问题定义及题目分析
1.求解线性方程是实际中常遇到的问题。
而一般求解方法有两种,一个是直接求解法,一个是迭代法。
2.直接法是就是通过有限步四则运算求的方程准确解的方法。
但实际计算中必然存在舍入误差,因此这种方法只能得到近似解。
3.迭代法是先给一个解的初始近似值,然后按一定的法则求出更准确的解,即是用某种极限过程逐步逼近准确解的方法。
4.这次实践,将采用直接法:列主元高斯消去法,全主元高斯消去法;迭代法:雅克比迭代法和高斯-赛德尔迭代法。
三 详细设计 设有线性方程组
11112211n n a x a x a x b +++=
21122222n n a x a x a x b +++= 1122n n nn n n a x a x a x b +++=
令A= 1112121
2221
2
n n n n nn a a a a
a a a a a ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥
⎣⎦ x=
12n x x x ⎡⎤⎢⎥
⎢⎥⎢⎥⎢⎥⎣⎦ b= 12n b b b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦
一. 上三角方程组的求解很简单。
所以若能将方程组化为上三角方程组,那就很容易得到方程的解,高斯消去法则是一种简洁高效的方法,可以将方程组化为上三角方程组,但是这种方法必须满足每一步时0kk a =才能求解,还有当kk a 绝对值很小时,作分母会引起较大的舍入误差。
因此在消元过程中应尽量选择绝对值较大的系数作为主元素。
1.列主元消去法很好的解决这一问题。
将1x 的系数1(1)k a k n ≤≤中绝对值最大者作为主元素,交换第一行和此元素所在的行,然后主元消去非主元项1x 的系数,。
一般的,在k x 的一些相关系数之前,先从,1,,,k k k k n k a a a + 中选出绝对值最大者作为主元素,交换第k 行和次主元素所在的行,再消去k x 的一些相关系数,只要方程组A 非奇异,则每个主元定不为0,故消元过程一定能将方程组化为上三角方程组。
将方程用增广矩阵[](1)()ij n n A b a ⨯+= 表示。
(1) 消元过程
对01,2,1k n =- ,;
1. 选主元,找{},1,k i k k n ∈+ ,使得,max ()k
i k ik a a k i n =≤≤
2. ,=0k
i k a ,则矩阵A 奇异,程序结束,否则执行下一步。
3. 如果k i k ≠,则交换第k 行与第k i 行对应的元素位置,即 (,1,,1)k
kj i j a a j k k n ↔=++
4. 消元,对1,2,,i k k n =++ 计算/ik ik kk l a a = 对1,2,,1j k k n =+++ 计算ij ij ik kj a a l a =- (2) 回代过程
,1/n n n nn x a a +=,对1,,2,1i n =- 计算 ,11()/n
i i n ij i ii j i x a a x a +=+=-∑
2. 全主元消去法(对列主元消去法的改进)
首先,从方程组系数ij a 中选出绝对值最大者作为主元,然后通过行交换及列交换,是主元素移到第一行第一列,然后消去第一列的非主元项的系数。
一般的,在第k 次消元之前,先从n-k+1行和n-k+1列中选取绝对值最大者作为主元素,交换第k 行和此主元素所在之行,再交换第k 列和此主元素所在之列。
然后消元,消去第k 列的一些系数,但每一步都必须跟踪记录12,n x x x 之间顺序的变化。
for(l=i;l<n;l++) 从剩余的n-i 列和n-i 列找主元 {
for(j=i;j<n;j++) 并将最大者记为A[k1][k2] if(fabs(A[l][j])>fabs(A[k1][k2]))k1=l,k2=j; }
if(k1!=i) 如果k1不等于i ,则第i 行与第k1行互换 for(j=0;j<=n;j++) {
p=A[i][j];
A[i][j]=A[k1][j]; A[k1][j]=p; }
for(i=0;i<n;i++)Y[i]=i; 初始化 Y[n]
if(k2!=i) if 成立则第i 列与第k2列互换 {
p=Y[i]; 相应的将Y[n]的第i 列与第k2列互换
Y[i]=Y[k2]; 则Y[n]起到了跟踪记录列交换的作用
Y[k2]=p;
for(j=0;j<n;j++) 第i 列与第k2列互换 {
p=A[j][i];
A[j][i]=A[j][k2]; A[j][k2]=p; } }
二 迭代法经常用于求解阶数高且矩阵稀疏的方程组。
1. 雅克比迭代法
设方程组系数矩阵A 的对角元素0ii a ≠(1,2,,)i n = ,max 为迭代的最大次数,e 为误差线。
(1)区初始向量000
12
(,,)n x x x x = ,令k=0; (2) 对1,2,,i n = ,计算
1
1,()/n
k k i
i ij j ii j j i
x
b a x a +=≠=-
∑
(3) 11
n
k k i i i x x e +=-<∑,则输出1k i x +,程序结束,否则进行下一步。
(4) 如果max k ≥,则并不收敛,程序结束,否则执行第2步。
2. 高斯-赛德尔迭代法(对雅克比迭代法的改进)
这种迭代法在于充分利用每一小步所求出的新分量。
因此在新分量一旦求出之时,马上、就用新分量代替。
即在雅克比迭代法法中求1
k i x +时用111121,k k k i x x x +++- 分别代替121,k k k
i x x x - 。
使得程序效率更高。
for(i=0;i<n;i++) for 循环进行迭代
{ 并记录每次迭代的误差 p=0;
for(j=0;j<n;j++) {
if(i==j) continue; 如果i 等于j 则不执行下一步 p=p+A[i][j]*X[j];
}
Y[i]=(A[i][n]-p)/A[i][i]; 新分量的计算
error=error+fabs(Y[i]-X[i]); 误差计算
X[i]=Y[i]; 新分量代替旧分量
}
四调试与测试分析
由于for循环嵌套语之多句及数组下标控制之难,在编程过程中困难重重。
因此就先将程序的大致算法写出来,再按照算法编程。
而C 语言数组下标从0 开始,因此按算法中所有下标减一编程。
在列主元消去法中,每次交换第i行和主元所在的第k行(按下标从1开始),没有交换第1,2…i-1列,但不影响结果。
而在编写全主元消去法时,每次能正确跟踪,结果总不对。
经过几次详细检查发现在交换行或列时,都是从第i行或第i列开始交换,发现列交换不能类比列主元消去法中的交换方法,必须考虑全局。
1.直接法测试用例:1 1 0 3 4
2 1 -1 1 1
3 -1 -1 3 -3
-1 2 3 -1 4
列主元高斯消去法输出全主元高斯消去法输出
X[1]=-1.33333 X[4]=1
X[2]=2.33333 X[3]=-0.333333
X[3]=-0.333333 X[1]=-1.33333
X[4]= 1 X[2]=2.33333
真实值为:X[1]=- 4/3,X[2]= 7/3,X[3]= -1/3, X[4]= 1。
2.迭代发测试用例:max=100
e=0.0001
10 -2 -1 3
-2 10 -1 15
-1 -2 5 10
雅克比迭代法输出高斯-赛德尔迭代法输出
迭代次数为12迭代次数为7
X[1]=0.999991 X[1]=0.999994
X[2]=1.99999 X[2]=2
X[3]=2.99999 X[3]=3
真实值为:X[1]=1,X[2]=2,X[3]=3。
两种直接法并无很大明显差别。
高斯-赛德尔迭代法明显优于雅克比迭代法。
五结果分析及总结
实验结果与预期一致,程序设计基本合理,能正确运行出结果。
在设计过程中能够及时发现问题并解决问题。
通过这次计算实践,对方程组求解的方法有了深入了解和掌握。
可视化实验结果为:
图1.列主元高斯消去法
图2.全主元高斯消去法
图3.雅克比迭代法
图4.高斯-赛德尔迭代法。