码迷,mamicode.com
首页 > 其他好文 > 详细

线性方程组的迭代解法——超松弛迭代法

时间:2019-12-30 14:34:32      阅读:103      评论:0      收藏:0      [点我收藏+]

标签:方程组   matlab   psi   ==   actor   radius   ros   func   break   

  1.代码

%%超松弛迭代法(此方法适用于大型稀疏矩阵但不适合与病态方程的解
%%线性方程组M*X = b,M是方阵,X0是初始解向量,epsilon是控制精度,omiga是松弛因子
function OIM = Overrelaxation_iterative_method(M,b,X0,epsilon)
[m,n] = size(M);
d = diag(M);L = zeros(m,n);U = zeros(m,n);D = zeros(m,n);eps = floor(abs(log(epsilon)));
ub = 10000;X = zeros(m,ucb);X(:,1) = X0;X_delta = X;X_end = zeros(m,1);k_end = 0;
for i = 1:m
    for j = 1:n
        if i > j
            L(i,j) = -M(i,j);
        elseif i < j
            U(i,j) = -M(i,j);
        elseif i == j
            D(i,j) = d(i);
        end
    end
end
J = D\(L+U);
MAX = Spectral_radius(J);
Optimum_relaxation_factor = 2/(1+sqrt(1-MAX^2));
disp(‘最佳松弛因子:‘);
abs(Optimum_relaxation_factor)
omiga = input(‘输入松弛因子:‘);
L_omiga = (D-omiga*L)\((1-omiga)*D+omiga*U);
f = omiga*((D-omiga*L)\b);
for k = 1:ub-1
    X(:,k+1) = L_omiga*X(:,k)+f;
    X_delta(:,k) = X(:,k+1)-X(:,k);
    delta = norm(X_delta(:,k),2);
    if delta < epsilon
        break
    end
end
disp(‘迭代解及迭代次数为:‘);
k
OIM = vpa([X(:,k)‘],eps);
%%谱半径
%%M是方阵
    function Sr = Spectral_radius(M)
        e = eig(M);
        Sr = max(e);
    end
end

  2.例子

clear all
clc
for i = 1:8
    for j = 1:8
        if i == j
            M(i,j) = 2.1;
        elseif i - j == 1
            M(i,j) = 1;
        elseif j - i == 1
            M(i,j) = -1;
        else
            M(i,j) = 0;
        end
    end
end
b = [1 2 3 4 4 3 2 1]‘;
X0 = [1 1 1 1 1 1 1 1]‘;
epsilon = 1e-4;

S = Overrelaxation_iterative_method(M,b,X0,epsilon)
M\b

  结果如下

最佳松弛因子:
ans =
    0.8540
输入松弛因子:0.8
迭代解及迭代次数为:
k =
    11
S =
[ 1.07163702, 1.25042379, 1.69753618, 1.81524709, 1.50955555, 0.985313468, 0.578713811, 0.200612427]
ans =
    1.0716
    1.2504
    1.6975
    1.8152
    1.5096
    0.9853
    0.5787
    0.2006

  

线性方程组的迭代解法——超松弛迭代法

标签:方程组   matlab   psi   ==   actor   radius   ros   func   break   

原文地址:https://www.cnblogs.com/guliangt/p/12119311.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!