matlab 求解线性方程组——迭代法

标签（空格分隔）： matlab 线性方程组 迭代 数值

文章目录

• • 标签（空格分隔）： matlab 线性方程组 迭代 数值
  • @[toc]
• matlab 求解线性方程组——迭代法
• • 1 测试迭代法求解函数GaussSeidelIteration()、JacobiIteration()
  • 2 迭代法求解函数GaussSeidelIteration()
  • 3 迭代法求解函数JacobiIteration()
  • 4 求解矩阵的范数matrixnorm()
  • 5 求解矩阵的范数条件数matrixcond()
  • 6 向量的p范数vectornorm()
  • 联系作者 definedone@163.com

matlab 求解线性方程组——迭代法

1 测试迭代法求解函数GaussSeidelIteration()、JacobiIteration()

clc,clear
% m=[1 2
%    -1 2]
% n1=matrixnorm([1 2;-1 2],1)
% n2=matrixnorm([1 2;-1 2],2)
% ninf=matrixnorm([1 2;-1 2],inf)
% 
% cond1=matrixcond(m,1)
% cond2=matrixcond(m,2)
% condinf=matrixcond(m,inf)

% a=[10.2 -0.25 -.3 
%     -.41 1.13 -.15 
%     -.25 -.14 1.21]
% b=[.515
%     1.555
%     2.78]
% precision=1e-3
% x=[2 2 2]'
% x=JacobiIteration(a,b,precision,x)
% xsystem=inv(a)*b

a=[6 -1 -1
    -1 6 -1
    -1 -1 6]
b=[11.33 32 42]'
precision=1e-3
x=[4 6 8]'
x=GaussSeidelIteration(a,b,precision,x) 
xsystem=inv(a)*b

运行结果：

iterationTimes =
     7
x =
    4.6661
    7.6189
    9.0475
xsystem =
    4.6661
    7.6189
    9.0475

2 迭代法求解函数GaussSeidelIteration()

function x=GaussSeidelIteration(a,b,precision,x) 
%时间：2013/12/08 姓名：邓能财   
disp('%GaussSeidel——迭代法求解线性方程组')
if nargin<4,  x=zeros(dim,1);
elseif nargin<3,  precision=1e-9; end
[dim,dim2]=size(a);
disp('迭代的方程为：x(k+1)=s*xk+f')
%错误的矩阵：
assert( dim==dim2 && dim==length(b),...
    ['Argument input error:   ',...
    'Matrix dimensions must agree.'])
d=diag(diag(a))
a=a-d;
l=-tril(a)
u=-triu(a)
invd_l=inv(d-l)
s=invd_l*u
f=invd_l*b
x_=ones(dim,1);
precision=precision/10;
x
iterationTimes=0
while any(abs(x-x_)>=precision)
    x_=x;
    x=s*x_+f
    iterationTimes=iterationTimes+1
end
end

3 迭代法求解函数JacobiIteration()

function x=JacobiIteration(a,b,precision,x)
%Jacobi迭代法求解线性方程组
%时间：2013/12/08 姓名：邓能财   
[dim,dim2]=size(a);
if nargin<4,  x=zeros(dim,1);
elseif nargin<3,  precision=1e-10; end
disp('迭代的方程为：x(k+1)=g*xk+f')
%错误的矩阵：
assert( dim==dim2 && dim==length(b),...
    ['Argument input error:   ',...
    'Matrix dimensions must agree.'])
d=diag(diag(a))
invd=inv(d)
g=invd*(d-a)
f=invd*b
x_=ones(dim,1);
iterationTimes=0
while any(abs(x-x_)>=precision)
    x_=x;
    x=g*x_+f
    iterationTimes=iterationTimes+1
end
end

4 求解矩阵的范数matrixnorm()

function n=matrixnorm(m,p) 
%求矩阵的p(1,2,inf)范数
%时间：2013/12/08 姓名：邓能财   
if nargin<2,  p=1; end
[dim,dim2]=size(m);
%错误的矩阵：
assert( dim==dim2,...
    ['Argument input error:   ',...
       'Matrix dimensions must agree.'])
switch p
    case 1
       disp('矩阵的一范数')
       sum_=sum(m,1)
       ssum_=sym(sum_)
       [max_,maxindex]=max(sum_)
       n=max_
    case inf
        disp('矩阵的无穷范数')
       sum_=sum(m,2)
       ssum_=sym(sum_)
       [max_,maxindex]=max(sum_)
       n=max_
     case 2
        disp('矩阵的2范数')
        t=m'*m
        st=sym(m'*m)
        %t的特征多项式
        syms x
        tx=x*eye(dim)-t
        fx=det(tx)
        eigenvalue=eig(t)
        seig=eig(st)
        sym(eigenvalue)
        real_eigenvalue=real(abs(eigenvalue))
        [max_,maxindex]=max(real_eigenvalue)
        n=sqrt(max_)
    otherwise
        error('第三个参数只能为1、2或inf')
end
end

5 求解矩阵的范数条件数matrixcond()

function c=matrixcond(m,p)
%时间：2013/12/08 姓名：邓能财   
if nargin<2,  p=1; end
[dim,dim2]=size(m);
%错误的矩阵：
assert( dim==dim2,...
    ['Argument input error:   ',...
    'Matrix dimensions must agree.'])
switch p
    case 1
        disp('矩阵的一范数条件数')
        inv_m=inv(m)
        sinv_m=inv(sym(m))
        c=matrixnorm(m,p)*matrixnorm(inv_m,p)
    case inf
        disp('矩阵的无穷范数条件数')
        inv_m=inv(m)
        sinv_m=inv(sym(m))
        c=matrixnorm(m,p)*matrixnorm(inv_m,p)
    case 2
        disp('矩阵的2范数条件数')
        inv_m=inv(m)
        sinv_m=inv(sym(m))
        c=matrixnorm(m,p)*matrixnorm(inv_m,p)
    otherwise
        error('第二个参数只能为1、2或inf')
end
end

6 向量的p范数vectornorm()

function n=vectornorm(x,p) %向量的p范数
%时间：2013/12/08 姓名：邓能财   
if nargin<2,  p=1; end
%错误的p：
assert( p>0,...
    ['Argument input error:   ',...
       'The next argument must be positive number! '])
 n=sum(abs(x).^p)^(1/p);
end

联系作者 definedone@163.com

end