一、数学原理
线性方程组 Ax=b 有惟一解,经过变换构造出一个等价同解方程组将上式改写成迭代式: , 选定初始向量 ,反复不断地使用迭代式,逐步逼近方程组的精确解,直到满足精度要求为止,这种方法称为迭代法。
对于给定的方程组可以构造各种迭代公式,但并非全部收敛 。其适用范围是迭代矩阵的谱半径小于1,即:。常用的迭代法有两种:雅克比迭代和Gauss–Saidel迭代。
雅克比迭代:
二、实验内容
分别计算Jacobi迭代法和Gauss–Saidel迭代法其迭代矩阵的谱半径,再用Gauss–Saidel迭代法求解下列线性方程组(3-1),用误差估计e≤10−6和最大迭代次数Nmax=500来控制迭代次数。
代码部分:
format compact
clc,clear;
e=1e-6;
Nmax=500;
A0=[4,3,0;3,4,-1;0,-1,4];
n=length(b0);
%Jacobi迭代法谱半径的计算
disp('Jacobi迭代法谱半径:')
D=diag(diag(A0));
pJ=max(abs(eig(inv(D)*(D-A0)))); )
fprintf('ρ(J)=%f\n',pJ);
%Gauss--Saidel迭代法谱半径的计算
disp('Gauss--Saidel迭代法谱半径:')
A1=tril(A0);
A2=triu(A0);
L=D-A1;
U=D-A2;
G=(inv(D-L))*U;
pG=max(abs(eig(G)));
fprintf('ρ(G)=%f\n',pG);
%Gauss--Saidel迭代法求解线性方程组,矩阵形式
x0=[0,0,0]'; %x0赋初值为(0,0,0)'
x=x0;x0=x+2*e;k=0;
b=b0';
A=A0;
A1=tril(A0);
iA1=inv(A1);
disp('Gauss--Saidel迭代法:')
while norm(x0-x,inf)>e&&k<Nmax
k=k+1;x0=x;
x=iA1*(A1-A0)*x0+iA1*b;
end
if k==Nmax
warning('已达到迭代上限次数,Gauss--Saidel迭代法失效');
else
x=x';
fprintf(' 解向量 迭代次数\n');
fprintf('%3.2f %3.2f %3.2f %d\n',x,k);
end
二、实验结果