【数值分析】线性方程组的迭代方法,jacobi,高斯赛德尔GS,SOR

线性方程组的迭代解法

2024年1月1日
#analysis



基本迭代法

Jacobi迭代

A = D − L − U A=D-L-U A=DLU
D x ( k + 1 ) = ( L + U ) x ( k ) + b Dx^{(k+1)}=(L+U)x^{(k)}+b Dx(k+1)=(L+U)x(k)+b
B j = D − 1 ( L + U ) = I − D − 1 A B_j =D^{-1}(L+U)=I-D^{-1}A Bj=D1(L+U)=ID1A
matlab实现

%% 迭代法例子
A = [2 -1 0;-1 3 -1;0 -1 2];
b = [1 8 -5]';
[x,i] = jacobi(A,b,1e-5,10000)
            
%% Jacobi迭代
% 输入矩阵A,向量b,精度,最大迭代次数
% 输出解向量x,迭代次数i
function [x,i] = jacobi(A,b,eps,max_iter)
    D = diag(diag(A));
    L = D-tril(A);
    U = D-triu(A);
    x = zeros(size(b)); %!
    for i = 1:max_iter
        x = D\(b+L*x+U*x);
        err = norm(b-A*x)/norm(b); %!
        if err<eps
            break;
        end
    end
end
高斯-赛德尔(GS)迭代

A = D − L − U A=D-L-U A=DLU
( D − L ) x ( k + 1 ) = U x ( k ) + b (D-L)x^{(k+1)}=Ux^{(k)}+b (DL)x(k+1)=Ux(k)+b
B g s = ( D − L ) − 1 U = I − ( D − L ) − 1 A B_{gs} =(D-L)^{-1}U=I-(D-L)^{-1}A Bgs=(DL)1U=I(DL)1A
matlab实现

%% 迭代法例子
A = [2 -1 0;-1 3 -1;0 -1 2];
b = [1 8 -5]';
[x,i] = GS(A,b,1e-5,10000)
            
%% GS迭代
% 输入矩阵A,向量b,精度,最大迭代次数
% 输出解向量x,迭代次数i
function [x,i] = GS(A,b,eps,max_iter)
    D = diag(diag(A));
    L = D-tril(A);
    U = D-triu(A);
    x = zeros(size(b)); %!
    for i = 1:max_iter
        x = (D-L)\(b+U*x);
        err = norm(b-A*x)/norm(b); %!
        if err<eps
            break;
        end
    end
end
SOR迭代

A = D − L − U A=D-L-U A=DLU
x ( k + 1 ) = x ( k ) + ω D − 1 ( L x ( k + 1 ) + U x ( k ) − D x ( k ) + b ) x^{(k+1)}=x^{(k)}+ \omega D^{-1}(Lx^{(k+1)}+Ux^{(k)}-Dx^{(k)}+b) x(k+1)=x(k)+ωD1(Lx(k+1)+Ux(k)Dx(k)+b)
x ( k + 1 ) = ( D − ω L ) − 1 [ ( 1 − ω ) D + ω U ] x ( k ) + ω ( D − ω L ) − 1 b x^{(k+1)}= (D- \omega L)^{-1}[(1- \omega )D+ \omega U]x^{(k)} + \omega (D- \omega L)^{-1}b x(k+1)=(DωL)1[(1ω)D+ωU]x(k)+ω(DωL)1b
B S O R = ( D − ω L ) − 1 [ ( 1 − ω ) D + ω U ] B_{SOR} = (D- \omega L)^{-1}[(1- \omega )D+ \omega U] BSOR=(DωL)1[(1ω)D+ωU]
matlab实现

%% 迭代法例子
A = [2 -1 0;-1 3 -1;0 -1 2];
b = [1 8 -5]';
[x,i] = SOR(A,b,1e-5,10000,1.1)
            
%% SOR迭代
% 输入矩阵A,向量b,精度,最大迭代次数
% 输出解向量x,迭代次数i
function [x,i] = SOR(A,b,eps,max_iter,w)
    D = diag(diag(A));
    L = D-tril(A);
    U = D-triu(A);
    x = zeros(size(b)); %!
    for i = 1:max_iter
        x = (D-w*L)\(((1-w)*D+w*U)*x + w*b);
        err = norm(b-A*x)/norm(b); %!
        if err<eps
            break;
        end
    end
end

迭代的收敛性分析和误差估计

排列矩阵 每行每列仅有唯一非零元的方阵。
可约矩阵 A {A} A n {n} n 阶矩阵, n ≥ 2 {n\ge2} n2 ,如果存在 n {n} n 阶排列矩阵 P {P} P ,使得
P T A P = [ A 11 A 12 0 A 22 ] P^ \mathrm TAP= \begin{bmatrix} A_{11} & A_{12} \\ 0 & A_{22} \end{bmatrix} PTAP=[A110A12A22]
其中 A 11 {A_{11}} A11 A 22 {A_{22}} A22 分别为 r {r} r 阶和 n − r {n-r} nr 阶方阵, 1 ≤ r ≤ n − 1 {1\le r\le n-1} 1rn1 ,则称 A {A} A 为可约矩阵,否则为不可约矩阵。
对角占优矩阵 A {A} A n {n} n 阶矩阵,满足
∣ a i i ∣ ≥ ∑ j = 1 , j ≠ i n ∣ a i j ∣    ,    i = 1 , 2 , ⋯   , n | a_{ii} |\ge \sum_{j=1,j\ne i}^{ n}|a_{ij}| \,\,,\,\, i=1,2,\cdots,n aiij=1,j=inaij,i=1,2,,n
即对角元素大于等于该行其他元素的和,如果 A {A} A 中至少有一行使不等式严格成立,则称A为弱对角占优矩阵,如果每一行都使不等式严格成立,则称 A {A} A 为严格行对角占优矩阵。

一些定理

  • 如果 n {n} n 阶矩阵 A {A} A 是严格对角占优矩阵或不可约弱对角占优矩阵,则 A {A} A 是非奇异矩阵
  • n {n} n 阶矩阵 A {A} A k {k} k 次幂 A k → 0 {A^k\to0} Ak0 的充要条件为谱半径 ρ ( A ) < 1 {\rho (A)<1} ρ(A)<1
  • 任一矩阵 A {A} A 的谱半径均不大于 A {A } A 的任一与某一向量范数相容的矩阵范数,即 ρ ( A ) ≤ ∣ ∣ A ∣ ∣ {\rho(A)\le ||A||} ρ(A)∣∣A∣∣
  • 对于迭代格式
    x ( k + 1 ) = B x ( k ) + g x^{(k+1)}=Bx^{(k)}+g x(k+1)=Bx(k)+g
    给定任意的初值 x ( 0 ) {x^{(0)}} x(0) ,有下列收敛结果和误差估计0
    1. 迭代格式收敛的充要条件为谱半径 ρ ( B ) < 1 {\rho(B)<1} ρ(B)<1
    2. 如果 ∣ ∣ B ∣ ∣ < 1 {||B||<1} ∣∣B∣∣<1 ,则有估计
      ∣ ∣ x ( k ) − x ∗ ∣ ∣ ≤ ∣ ∣ B ∣ ∣ k 1 − ∣ ∣ B ∣ ∣ ∣ ∣ x ( 1 ) − x ( 0 ) ∣ ∣ ∣ ∣ x ( k ) − x ∗ ∣ ∣ ≤ ∣ ∣ B ∣ ∣ 1 − ∣ ∣ B ∣ ∣ ∣ ∣ x ( k ) − x ( k − 1 ) ∣ ∣ \begin{align*} ||x^{(k)}-x ^{*} ||\le& \frac{||B||^k}{1-||B||}||x^{(1)}-x^{(0)}|| \\ \\ ||x^{(k)}-x ^{*} ||\le& \frac{||B||}{1-||B||}||x^{(k)}-x^{(k-1)}|| \end{align*} ∣∣x(k)x∣∣∣∣x(k)x∣∣1∣∣B∣∣∣∣Bk∣∣x(1)x(0)∣∣1∣∣B∣∣∣∣B∣∣∣∣x(k)x(k1)∣∣
  • A {A} A 是严格对角占优矩阵或不可约弱对角占优矩阵,则Jacobi迭代和GS迭代都收敛
  • A {A} A 对称正定,则Jacobi迭代收敛的充要条件为 2 D − A {2D-A} 2DA 也是对称正定矩阵
  • SOR迭代收敛的必要条件为 0 < ω < 2 {0< \omega <2} 0<ω<2
  • 系数矩阵 A {A} A 对称正定,则 0 < ω < 2 {0<\omega <2} 0<ω<2 时SOR迭代收敛

例题看同济《现代数值计算》习题6.6。


下链


  • 23
    点赞
  • 32
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
首先,我们需要明确雅可比迭代法和高斯德尔迭代法的迭代公式: 雅可比迭代法: $$ x_i^{(k+1)}=\frac{1}{a_{ii}} \left(b_i-\sum_{j=1,j\neq i}^{n} a_{ij}x_j^{(k)}\right),\quad i=1,2,\cdots,n $$ 高斯德尔迭代法: $$ x_i^{(k+1)}=\frac{1}{a_{ii}} \left(b_i-\sum_{j=1}^{i-1} a_{ij}x_j^{(k+1)}-\sum_{j=i+1}^{n} a_{ij}x_j^{(k)}\right),\quad i=1,2,\cdots,n $$ 其中,$a_{ij}$ 和 $b_i$ 分别表示系数矩阵和常数向量的元素,$x_i^{(k)}$ 表示第 $k$ 次迭代后第 $i$ 个未知量的值。 下面是 MATLAB 的实现: % 雅可比迭代法 function [x, k] = Jacobi(A, b, x0, epsilon, maxIter) % A:系数矩阵 % b:常数向量 % x0:迭代初值 % epsilon:迭代精度 % maxIter:最大迭代次数 % x:方程组的解 % k:迭代次数 n = length(b); x = x0; for k = 1 : maxIter for i = 1 : n temp = b(i); for j = 1 : n if j ~= i temp = temp - A(i, j) * x(j); end end x(i) = temp / A(i, i); end if norm(A * x - b, inf) < epsilon return; end end % 高斯德尔迭代法 function [x, k] = GaussSeidel(A, b, x0, epsilon, maxIter) % A:系数矩阵 % b:常数向量 % x0:迭代初值 % epsilon:迭代精度 % maxIter:最大迭代次数 % x:方程组的解 % k:迭代次数 n = length(b); x = x0; for k = 1 : maxIter for i = 1 : n temp = b(i); for j = 1 : i - 1 temp = temp - A(i, j) * x(j); end for j = i + 1 : n temp = temp - A(i, j) * x(j); end x(i) = temp / A(i, i); end if norm(A * x - b, inf) < epsilon return; end end

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值