线性方程组的数值解法
线性方程组数值方法
线性方程组数值方法Ax ，必须满足系数矩阵的对角元素非零的条件。
1、回代：用回代法求解上三角线性方程组bfunction x=backsub(A,b)%Input -A is an nxn upper-triangular nonsingular matrix% -b is an nx1 matrix%Output-x is the solution to the linear system Ax=b%Find the dimension of b and initialize xn=length(b);x=zeros(n,1);x(n)=b(n)/A(n,n);for k=n-1:-1:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);end例：A=[1 2 1 4;0 -4 2 -5;0 0 -5 -7.5;0 0 0 -9];b=[13;2;-35;-18];x=backsub(A,b)2、高斯消去法(不选主元素)function x=gauss(A,b)%Input -A is an nxn nonsingular matrix% -b is an nx1 matrix%Output-x is an nx1 matrix containing the solution to Ax=b%Initialize x[n,n]=size(A);x=zeros(n,1);%Form the augmented matrix:Aug=[A|b]Aug=[A b];for k=1:n-1%Elimination process for column kfor i=k+1:nm=Aug(i,k)/Aug(k,k);Aug(i,k:n+1)=Aug(i,k:n+1)-m*Aug(k,k:n+1);endend%Output Elimination matrixAug%Back SubstitutionA=Aug(1:n,1:n);b=Aug(1:n,n+1);x(n)=b(n)/A(n,n);for k=n-1:-1:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);end例： A=[4 -9 2;2 -4 6;-1 1 -3];b=[5;3;-4]; x=gauss(A,b)3、列主元消去法function x=gauss_column(A,b)%Input -A is an nxn nonsingular matrix% -b is an nx1 matrix%Output-x is an nx1 matrix containing the solution to Ax=b%Initialize x and the temporary storage matrix C[n,n]=size(A);x=zeros(n,1);C=zeros(1,n+1);%Form the augmented matrix:Aug=[A|b]Aug=[A b];for k=1:n-1%Partial pivoting for column k[y,j]=max(abs(Aug(k:n,k)));%Interchange row k and jC=Aug(k,:);Aug(k,:)=Aug(j+k-1,:);Aug(j+k-1,:)=C;if Aug(k,k)==0'A was singular,No unique solution'breakend%Elimination process for column kfor i=k+1:nm=Aug(i,k)/Aug(k,k);Aug(i,k:n+1)=Aug(i,k:n+1)-m*Aug(k,k:n+1);endend%Back Substitution on [U|y] using backsubx=backsub(Aug(1:n,1:n),Aug(1:n,n+1));例：A=[1 2 1 4;2 0 4 3;4 2 2 1;-3 1 3 2];b=[13;28;20;6];x=gauss_column(A,b)4、LU分解法(不选主元素)function x=lufact(A,b)%Input -A is an nxn matrix% -b is an nx1 matrix%Output-x is an nx1 matrix containing the solution to Ax=b%Initilize x,y,[n,n]=size(A);x=zeros(n,1);y=zeros(n,1);for k=1:n-1%Calculate multiplier and place in subdiagonal portion of A for i=k+1:nmult=A(i,k)/A(k,k);A(i,k)=mult;A(i,k+1:n)=A(i,k+1:n)-mult*A(k,k+1:n);endend%Output LU factorization compact matrixA%solve for yy(1)=b(1);for i=2:ny(i)=b(i)-A(i,1:i-1)*y(1:i-1);end%Solve for xx(n)=y(n)/A(n,n);for i=n-1:-1:1x(i)=(y(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);end例：A=[4 -9 2;2 -4 6;-1 1 -3];b=[5;3;-4];x=lufact(A,b)5、列主元LU分解法（PA＝LU）A是非奇异矩阵function x=lu_column(A,b)%Input -A is an nxn matrix% -b is an nx1 matrix%Output-x is an nx1 matrix containing the solution to Ax=b%Initilize x,y,the temporary storage matrix C,and the row%permutation information matrix R[n,n]=size(A);x=zeros(n,1);y=zeros(n,1);C=zeros(1,n);R=1:n;for k=1:n-1%Find the pivot row for column k[maxl,j]=max(abs(A(k:n,k)));%Interchange row k and jC=A(k,:);A(k,:)=A(j+k-1,:);A(j+k-1,:)=C;d=R(k);R(k)=R(j+k-1);R(j+k-1)=d;if A(k,k)==0'A is singular,No unique solution'breakend%Calculate multiplier and place in subdiagonal portion of A for i=k+1:nmult=A(i,k)/A(k,k);A(i,k)=mult;A(i,k+1:n)=A(i,k+1:n)-mult*A(k,k+1:n);endend%Output LU factorization compact matrixA%solve for yy(1)=b(R(1));for i=2:ny(i)=b(R(i))-A(i,1:i-1)*y(1:i-1);end%Solve for xx(n)=y(n)/A(n,n);for i=n-1:-1:1x(i)=(y(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);end例：A=[1 2 6;4 8 -1;-2 3 -5];b=[1;2;3];x=lu_column(A,b)6、Jacobi迭代:程序可用的充分条件是A具有严格对角优势。
线性代数方程组的数值解法讨论
线性代数方程组的数值解法讨论解线性方程组的方法，主要分为直接方法和迭代方法两种。
直接法是在没有舍入误差的假设下能在预定的运算次数内求得精确解。
而实际上，原始数据的误差和运算的舍入误差是不可以避免的，实际上获得的也是近似解。
迭代法是构造一定的递推格式，产生逼近精确解的序列。
对于高阶方程组，如一些偏微分方程数值求解中出现的方程组，采用直接法计算代价比较高，迭代法则简单又实用，因此比较受工程人员青睐。
小组成员本着工程应用，讨论将学习的理论知识转变为matlab 代码。
讨论的成果也以各种代码的形式在下面展现。
1 Jacobi 迭代法使用Jacobi 迭代法，首先必须给定初始值，其计算过程可以用以下步骤描述： 步骤1 输入系数矩阵A ，常熟向量b ，初值(0)x ，误差限ε，正整数N ，令1k =.步骤2 (0)11ni i ij jj ii j i x b a x a =≠⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦∑，(0)j x 代表(0)x 的第j 个分量。
步骤3 计算11ni i ij j j ii j i y b a x a =≠⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦∑，判断1max i i i n x y ε≤≤-<，如果是，则结束迭代，转入步骤5；否则，转入步骤4。
步骤4 判断k N =？如果是，则输出失败标志；否则，置1k k =+，i i x y ⇐，1,2,,i n =，转入步骤2。
步骤5 输出12,,n y y y 。
雅可比迭代代码function [x,k]=Fjacobi(A,b,x0,tol)% jacobi 迭代法 计算线性方程组% tol 为输入误差容限，x0为迭代初值max1= 300; %默认最多迭代300，超过要300次给出警告 D=diag(diag(A)); L=-tril(A,-1);U=-triu(A,1); B=D\(L+U); f=D\b; x=B*x0+f;k=1; %迭代次数while norm(x-x0)>=tol x0=x;x=B*x0+f; k=k+1;if(k>=max1)disp('迭代超过300次，方程组可能不收敛'); return; end%[k x'] %显示每一步迭代的结果 End2 高斯赛德尔迭代由Jacobi 迭代法中，每一次的迭代只用到前一次的迭代值，若每一次迭代充分利用当前最新的迭代值，即在计算第i 个分量(1)k i x +时，用最新分量11()k x +，12()k x +…(1)1k i x +-代替旧分量)1(k x ', )2(k x …)3(k x 就得到高斯赛德尔迭代格式，其数学表达式为：1(1)(1)()111(1,2,,)i n k k k ii ij j ij j j j i ii xb a x a x i n a -++==+⎛⎫=--= ⎪⎝⎭∑∑具体形式如下：()()()(1)()()()11221331111(1)(1)()()22112332222(1)(1)(1)(1)(1)112233,11111k k k k n n k k k k n n k k k k k n n n n n n n n nnx a x a x a x b a x a x a x a x b a x a x a x a x a x b a ++++++++--=----+=----+⋯⋯⋯⋯⋯⋯=-----+矩阵形式表示为：()(1)1(1)()(0,1,2,,),k k k k n +-+=++=x D Lx Ux b将(1)(1)()(0,1,2,,)k k k k n ++=++=Dx Lx Ux b 移项整理得： (1)1()1()()(0,1,2,,))k k x D L Ux D L b k n +--=-+-=记11(),()--=-=-M D L U g D L b ，则(1)()k k x x +=+M g高斯塞德尔迭代function [x,k]=Fgseid(A,b,x0,tol)%高斯-塞德尔迭代法 计算线性方程组 % tol 为误差容限max1= 300; %默认最高迭代300次D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); G=(D-L)\U; f=(D-L)\b; x=G*x0+f;k=1; while norm(x-x0)>=tol x0=x;x=G*x0+f; k=k+1;if(k>=max1)disp('迭代次数太多，可能不收敛'); return; end% [k,x'] %显示每一步迭代结果 End3 超松弛迭代法在工程中最常遇到的问题便是线性代数方程组的求解,而线性代数方程组的求解一般可以分为两类,一类是直接法(精确法),包括克莱姆法则方法、LD 分解法等，另一类是迭代法（近似法），包括雅克比迭代法、高斯迭代法、超松弛迭代法等。
数值分析计算方法实验报告
end;
end;
X=x;
disp('迭代结果:');
X
format short;
输出结果：
因为不收敛，故出现上述情况。
4.超松弛迭代法：
%SOR法求解实验1
%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]'
b=[10,5,-2,7]'
[m,n]=size(A);
if m~=n
error('矩阵A的行数和列数必须相同');
return;
end
if m~=size(b)
error('b的大小必须和A的行数或A的列数相同');
return;
end
if rank(A)~=rank([A,b])
error('A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');
3．实验环境及实验文件存档名
写出实验环境及实验文件存档名
4.实验结果及分析
输出计算结果，结果分析和小结等。
解：1.高斯列主元消去法：
%用高斯列主元消去法解实验1
%高斯列主元消元法求解线性方程组Ax=b
%A为输入矩阵系数,b为方程组右端系数
%方程组的解保存在x变量中
format long;
A=[2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13]
return;
end
c=n+1;
A(:,c)=b;
for k=1:n-1
高斯顺序消去法的条件 -回复
高斯顺序消去法的条件-回复高斯顺序消去法，也被称为高斯消元法或高斯消元算法，是一种用于解线性方程组的数值解法。
它的基本思想是通过一系列的行变换将线性方程组转化为一个上三角矩阵，然后通过回代求解出方程的解。
高斯顺序消去法适用于满足以下条件的线性方程组：1. 方程个数等于未知数个数：高斯顺序消去法只适用于方程个数等于未知数个数的情况。
假设方程个数为n，那么未知数个数也应为n。
2. 矩阵的行列式不为0：在进行高斯顺序消去法时，需要将系数矩阵进行行变换，其中包括将某一行乘以一个非零常数，或将某一行加到另一行上。
若矩阵的行列式为0，说明矩阵的行线性相关，无法进行行变换得到上三角矩阵。
3. 未知数的个数必须大于或等于方程组的秩：方程组的秩指的是矩阵的秩，即线性方程组的系数矩阵经过高斯消元法得到的上三角矩阵中非零行的个数。
未知数个数必须大于或等于方程组的秩，否则方程组存在自由变量，无法使用高斯顺序消去法求解。
4. 系数矩阵的系数不应过大或过小：在进行高斯顺序消去法时，系数矩阵中的系数大小对计算结果的精确性有影响。
若系数值过大或过小，会导致计算过程中出现大数吃小数或小数吃小数的情况，从而降低计算结果的准确性。
基于上述条件，以下是高斯顺序消去法的具体步骤：步骤1：构造增广矩阵将线性方程组的系数矩阵与常数项矩阵按列连接起来，形成一个增广矩阵。
步骤2：主元确定找到增广矩阵的第一列中绝对值最大的元素（称为主元），将其所在行与第一行进行交换。
步骤3：主元归一化将主元所在行的元素都乘以一个适当的常数，使主元变为1。
步骤4：消元过程从第一行的下一行开始，将各行的第一列元素消为零。
具体操作是，将下一行第一列的元素乘以一个适当的常数，然后将得到的结果加到第一行的对应的元素上。
重复这一步骤，直到第一列的所有元素都变为零。
步骤5：重新确定主元在步骤4的过程中可能会存在一些舍入误差，导致第一列中的某些元素不为零。
因此，在消元过程结束后，重新找到第二列中的主元，将其所在行与第二行进行交换，并将主元归一化。
数值分析-北交大-王兵团-3-线性方程组解法 (1)
©
追赶法求解公式为：
追赶法算法
用追赶法来求解三对角线性方程组， 计算量只是5n-4，这比Gauss消元法的计算 量要小很多。
©
第3章 线性方程组解法
§3.5 线性方程组解对系数的敏感性
©
1、解对系数敏感性的相对误差 设方程组Ax=b的解为
扰动方程组的准确解为
有
©
用上述过程求解 的方法称为追赶法解法。
©
定理3.7
Sor法收敛的必要条件是松弛因子满足0<<2 证明
©
2、误差估计 定理3.8 设矩阵B的某种矩阵范数
证明参照非线性方程求根定理的证明， 将：绝对值换成范数、函数换成矩阵，注意范数关系 的使用，
©
例3.1 用Jacobi 迭代法解线性方程组 解
Jacobi迭代收敛！
故所求近似解为 准确解：
©
第3章 线性方程组解法
§3.4 线性方程组的直接解法
©
一、Gauss消元法 1、基本思想 先将线性方程组通过消元方法化为同解的上三角
方程组，然后从该三角方程组中按第n个方程、第n1个方程、…、第1个方程的顺序，逐步回代求出线 性方程组的解。
2、构造原理 Gauss消元法的求解过程分为两个: “消元”：把原方程组化为上三角方程组； “回代”：求上三角方程组的解。
©
计算量
©
2）Gauss消元法矩阵解释 第1步消元
第n-1步消元后，有
©
L是下三角阵，U是上 三角阵。
A=D-L-U ?
例：研究线性方程组
的Gauss消元法求解结果，假设计算在4位浮点十进 制数的计算机上求解。
解：
用Gauss消元法得
©
用Gauss消元法求解得 其准确解为
数值计算方法第3章解线性方程组的数值解法1
,i
2 ,3 ,...,
n
a
(1 11
)
A（ 1） A ( 2 )
a (1) 11
a (2) 22
...... ......
......
a (2) n2
......
a a
(1) 1n
(2) 2n
a
(2 nn
)
b (1)
b (2)
[
b
( 1
1
)
b (2) 2
a(k) kk
...
a(k) kn
... ... ...
...
...
a(n) nn
b1(1) b2(2)
...
bk(k)
...
bn(n)
21
高斯顺序消去法
也就是对于方程组AX=b系数矩阵做：
ai(jkl1i)k
a(k) ik
a(k) ij
/
a(k) kk
3)顺序消元
31
高斯列主元消去法
第k步
从A ( k ) 的第
k
列
a (k) kk
，a (k) k 1k
，...a
(k) nk
中选取绝对值
最大项，记录所在行，即
|a(k) ikk
|m kina|axi(kk)
|
记 lik
若 l k 交换第k行与l行的所有对应元素，再 进行顺序消元。
32
其中， lii 0, i 1,2,..., n
（1）
10
高斯顺序消元法
线性代数方程组的数值解法_百度文库
线性代数方程组的数值解法【实验目的】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 线性方程组的数值解法1．1 题目的主要研究内容及预期达到的目标（1）实现雅可比迭代算法（2）求得给定线性方程组的解1．2 题目研究的工作基础或实验条件（1）笔记本（2）C-Free 5应用开发平台1．3 设计思想将线性方程组的系数矩阵A分成三部分D、L和U,其中满足A=D-L-U(具体关系参考教材P188)。
由Ax=b和雅可比迭代法思想可以推导得出最终的计算公式，然后根据公式和算法思想即可编写程序实现。
1．4 流程图1．5 主要程序代码(要求必须有注释)#include<iostream>#include<Windows.h>using namespace std;void JacobiIteration(){int n,k;cout<<"输入你将输入的系数矩阵阶数n=";cin>>n;double e,w[n];double Matrix[n][n],b[n];ZeroMemory(Matrix,sizeof(Matrix));//初始化系数矩阵，元素置零cout<<"输入系数矩阵(类似：\na1 a2 a3\nb1 b2 b3\nc1 c2 c3):"<<endl;for(int i=0;i<n;++i)for(int j=0;j<n;++j)cin>>Matrix[i][j];cout<<"依次输入"<<n<<"个线性方程组的右式值：";for(int i=0;i<n;++i)cin>>b[i];double x[n];for(int i=0;i<n;++i)x[i]=0;cout<<'\n'<<"输入迭代次数k=";cin>>k;for(int i=0;i<k;++i){for(int m=0;m<n;++m)w[m]=0;//清零，以便记录累加和for(int m=0;m<n;++m){for(int l=0;l<n;++l){if(l!=m)w[m]+=Matrix[m][l]*x[l];}}for(int j=0;j<n;++j){x[j]=(b[j]-w[j])/Matrix[j][j];}}cout<<"该线性方程组经雅可比迭代法求解得："<<endl;for(int i=0;i<n;++i){cout<<"x"<<i+1<<"="<<x[i]<<endl;//输出最终求解结果}}int main(){JacobiIteration();}1．6 运行结果及分析所用测试范例为教材P180的例1，最终求得结果同教材所给结果相同。
线性方程组数值解法总结
好久没来论坛，刚刚发现以前的帖子现在那么火很欣慰，谢谢大家支持！今天趁着不想做其他事情，把线性方程组的数值解法总结下，有不足的地方希望大神指教！数学建模中也会用到线性方程组的解法，你会发现上10个的方程手动解得话把你累个半死，而且不一定有结果，直接用matlab的函数，可以，关键是你不理解用着你安心吗？你怎么知道解得对不对？我打算开个长久帖子，直到讲完为止！这是第一讲，如有纰漏请多多直接，大家一起交流！线性方程组解法有两大类：直接法和迭代法直接法是解精确解，这里主要讲一下Gauss消去法，目前求解中小型线性方程组（阶数不超过1000），它是常用的方法，一般用于系数矩阵稠密，而有没有特殊结构的线性方程组。
首先，有三角形方程组的解法引入Gauss消去法，下三角方程组用前代法求解，这个很简单，就是通过第一个解第二个，然后一直这样直到解出最后一个未知数，代码如下：前代法：function [b]= qiandai_method(L,b)n=size(L,1); %n 矩阵L的行数for j=1:n-1 %前代法求解结果存放在b中b(j)=b(j)/L(j,j);b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j);endb(n)=b(n)/L(n,n);上三角方程组用回代法，和前面一样就是从下面开始解x，代码：后代法：function [y]=houdai_method(U,y)n=size(U,1); %n 矩阵L的行数for j=n:-1:2 %后代法求解结果存放在y中y(j)=y(j)/U(j,j);y(1:j-1)=y(1:j-1)-y(j)*U(1:j-1,j);endy(1)=y(1)/U(1,1);Gauss消去的前提就是这两个算法：具体思想是把任何一个线性方程组的系数矩阵A，分解为一个上三角和一个下三角的乘积，即A=LU，其中L为下三角，U为上三角。
那么具体怎么做呢？有高斯变换，什么是高斯变换？由于时间有限我不可能去输入公式，所以我用最平白的话把它描述出来。
数值分析-线性方程组的直接解法
算法 Gauss(A,a,b,n,x)
1. 消元 For k=1,2, … , n-1 1.1 if akk=0 , stop; 1.2 For i=k+1,k+2, …, n 1.2.1 l ik=aik /akk => aik 1.2.2 For j=k+1，k+2, … ,n ai j -aik ak j =>aij 1.2.3 bi -aik bk=> bi 2. 回代 2.1 bn / an=>xn; 2.2 For i=n-1,n-2, …, 2,1 2.2.1 bk => S 2.2.2 For j=k+1，k+2, … ,n S –akj xj =>S 2.2.3 S/ akk => xk a1 1 a1 2 a13 a2 1 a2 2 a23
线性方程组的直接解法
刘 斌
线性方程组的直接解法
§1 Gauss消去法 1.1 顺序Gauss消去法
1.2
§2 2.1 2.2 2.3
列主元Gauss消去法
Gauss消去法的矩阵运算 Doolittle分解法 平方根法
直接三角分解方法
2.4
追赶法
引入
在科学计算中，经常需要求解含有n个未知量 的n个方程构成的线性方程组 a11 x1 a12 x2 a1n xn b1 a21 x1 a22 x2 a2 n xn b2 (1) an1 x1 an 2 x2 ann xn bn
(1) a12 ( 2) a22 0
(1) (1) a13 a1 n ( 2) ( 2) a23 a2 n ( 3) ( 3) a33 a3 n
0
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
贵州师范大学数学与计算机科学学院学生实验报告
课程名称： 数值分析 班级： 实验日期：
学 号： 姓名： 指导教师：
实验成绩：
一、实验名称
实验五:线性方程组的数值解法
二、实验目的及要求
1. 让学生掌握用列主元gauss消去法、超松弛迭代法求解线性方程组.
2. 培养Matlab编程与上机调试能力.
三、实验环境
每人一台计算机，要求安装Windows XP操作系统，Microsoft office2003、
MATLAB6.5(或7.0).
四、实验内容
1. 编制逐次超松弛迭代(SOR迭代)函数(子程序),并用于求解方程组
141414144321432143214321xxxx
xxxx
xxxx
xxxx
取初始向量Tx)1,1,1,1()0(,迭代控制条件为 5)1()(1021||||kkxx
请绘制出迭代次数与松弛因子关系的函数曲线,给出最佳松弛因子.SOR迭代
的收敛速度是否一定比Gauss-Seidel迭代快?
2. 编制列主元 Gauss 消去法函数(子程序),并用于解
615318153312321321321xxx
xxx
xxx
要求输出方程组的解和消元后的增广矩阵.
注：题2必须写实验报告
五、算法描述及实验步骤
第一步：
先编出列主元高斯消去法的程序
Function [Ab,x]=liegauss(Ab,n)
k=1;
while k
if abs(Ab(i,k))>abs(Ab(k,k))
b=Ab(i,:);Ab(i,:)=Ab(k,:);Ab(k,:)=b;
end
end
for i=k+1:n
Ab(i,k:n+1)=-Ab(i,k)/Ab(k,k)*Ab(k,k:n+1)+Ab(i,k:n+1);
end
k=k+1;
end
x=zeros(n,1);
x(n)=Ab(n,n+1)/Ab(n,n);
for i=n-1:-1:1
x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);
end
第二步：
在MATLAB命令窗口里输入上述方程组的增广矩阵
Ab=[12,-3,3,15;-18,3,-1,-15;1,1,1,6]，输入n的值，调用所编的程序。
六、调试过程及实验结果
>> Ab=[12,-3,3,15;-18,3,-1,-15;1,1,1,6];
>> n=3;
>> [Ab,x]=liegauss(Ab,n)
Ab =
-18.0000 3.0000 -1.0000 -15.0000
0 1.1667 0.9444 5.1667
0 0 3.1429 9.4286
x =
1.0000
2.0000
3.0000
>>
七、总结
列主元消去法是数值稳定的方法，计算量小，可以避免绝对值小的数作除数，从
而避免舍入误差的恶性增长，所以实际计算中主要采用列主元消去法来解决中小
规模的线性方程组和某些大型稀疏线性方程组。
八、附录（源程序清单）
Function [Ab,x]=liegauss(Ab,n)
k=1;
while k
if abs(Ab(i,k))>abs(Ab(k,k))
b=Ab(i,:);Ab(i,:)=Ab(k,:);Ab(k,:)=b;
end
end
for i=k+1:n
Ab(i,k:n+1)=-Ab(i,k)/Ab(k,k)*Ab(k,k:n+1)+Ab(i,k:n+1);
end
k=k+1;
end
x=zeros(n,1);
x(n)=Ab(n,n+1)/Ab(n,n);
for i=n-1:-1:1
x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);
end