MATLAB应用求线性方程组的通解

理解线性方程组直接法与迭代法思想，掌握常用算法的设计，掌握用matlab实现的数值解法。

1、编写列主元消去法程序，并举例子。编写ＬＵ分解法程序，并举例子。对两种算法作出对比。利用MATLAB函数“”计算上面例子。

2、编写Jacobi或Gauss-Seidel迭代法程序，并举例子。编写SOR迭代法程序，并举例子选择较佳的松弛因子。对比说明两种迭代法的收敛速度。

3、谈一谈你对线性方程组直接法与迭代法的理解。

1、利用列主元法解方程：
gaussh.m文件为：

function x=gauss(a,b)

%%判断是否可以求出解

temp1=size(a);temp2=size(b);

if temp1(1)~=temp1(2)|temp1(1)~=temp2(1)|temp2(2)~=1|det(a)==0

    error('矩阵或向量的大小不对应或矩阵为奇异矩阵')

    return;

end

%%化为三角矩阵过程

leng=temp1(1);

for i=1:leng-1

     pos=i;num=abs(a(i,i));

    for j=i+1:leng

        if abs(a(j,i))>num

            num=abs(a(j,i));

            pos=j;

        end

    end

    if pos~=i

        for k=i:leng

            temp=a(i,k);

            a(i,k)=a(pos,k);

            a(pos,k)=temp;

        end

        temp=b(pos);

        b(pos)=b(i);

        b(i)=temp;

    end

       for j=i+1:leng

        temp=a(j,i)/a(i,i);

        for h=i:leng

            a(j,h)=a(j,h)-a(i,h)*temp;

        end

        b(j)=b(j)-temp*b(i);

    end

end

%%进行回代

x(leng)=b(leng)/a(leng,leng);

if leng>1

   for i=leng-1:-1:1

       x(i)=b(i);

       for j=i+1:leng;

           x(i)=x(i)-a(i,j)*x(j);

       end

       x(i)=x(i)/a(i,i);

   end

end

    x=x';

控制命令为：

>> clear

>> a=[1 2 3;2 5 2;3 1 5];b=[14 18 20]';

>> gaussh(a,b)

ans =
    1.0000

    2.0000

    3.0000

利用LU列主元法解方程
mylu.m文件内容为：

function [l,u,p]=mylu(a)

siz=size(a);

if siz(1)~=siz(2)

    error('矩阵非为方阵')

    return
end

if det(a)==0
    error('矩阵不能被三角分解')
end

si=siz(1);
u=a;p=eye(si);l=eye(si);

for i=1:si
    for j=i:si
        t(j)=u(j,i);
        for k=1:i-1
            t(j)=t(j)-u(j,k)*u(k,i);
        end
    end
    pos=i;ma=abs(t(i));
    for j=i+1:si
        if ma             ma=abs(t(j));

            pos=j;
        end
    end

    if pos~=i

        for j=1:si

            temp=u(i,j);

            u(i,j)=u(pos,j);

            u(pos,j)=temp;

        end

         for j=1:si

            temp=p(i,j);

            p(i,j)=p(pos,j);

            p(pos,j)=temp;

         end

       temp=t(pos);

       t(pos)=t(i);

       t(i)=temp;

    end

    u(i,i)=t(i);

    for j=i+1:si

        u(j,i)=t(j)/t(i);

    end

    for j=i+1:si

        for k=1:i-1

            u(i,j)=u(i,j)-u(i,k)*u(k,j);

        end

    end

end

l=tril(u,-1)+eye(si);

u=triu(u,0);

控制命令为：

a=[1 2 3;2 5 2;3 1 5];b=[14 18 20]';

>> [l,u,p]=mylu(a)

l =
    1.0000         0         0

    0.6667    1.0000         0

    0.3333    0.3846    1.0000
u =
        3.0000    1.0000    5.0000

         0    4.3333   -1.3333

         0         0    1.8462

p =
     0     0     1

     0     1     0

     1     0     0

>> p*b
ans =
    20

    18

    14
y=lans

y =
   20.0000

    4.6667

    5.5385
>> x=uy
x =
    1.0000

    2.0000

    3.0000

利用‘’直接求解为：

> x=ab
x =
    1.0000

    2.0000

    3.0000

2、用Jacobi解方程：

…………①

Jacobi.m文件为：

function [x,n]=jacobi(a,b,E)

temp1=size(a);temp2=size(b);

if temp1(1)~=temp1(2)|temp1(1)~=temp2(1)|temp2(2)~=1|det(a)==0

    error('矩阵或向量的大小不对应或矩阵为奇异矩阵')

    return;

end

%%化为三角矩阵过程

leng=temp1(1);N=-tril(a,-1)-triu(a,1);M=a+N;

t=1;

for i=1:leng

    t=t*M(i,i);

end

if t==0

    error('对角元素都为零！')

end

x=zeros(leng,1);

t=0;x=zeros(leng,1);

x2=N*x+b;

for i=1:leng  

    x1=x2(i)/M(i,i);

    t=t+(x1-x(i))^2;

    x(i)=x1;
end

e=E^2;

n=1;

while t>e

   t=0;

   x2=N*x+b;

   for i=1:leng     

       x1=x2(i)/M(i,i);

       t=t+(x1-x(i))^2;

       x(i)=x1;

   end

   n=n+1;

   if t>1000

       break;

       error('不收敛')

   end

end

控制窗口命令为：

a=[8 -3 2;4 11 -1; 6 3 12];b=[20 33 36]';

>> [x,n]=jacobi(a,b,0.001)

x =
    3.0003

    1.9999

    0.9997

n =
     9

用Gauss-Seidel迭代法解方程①

Gauss_seidel.m文件内容为：

function [x,n]=gauss_seidel(a,b,E)

temp1=size(a);temp2=size(b);

if temp1(1)~=temp1(2)|temp1(1)~=temp2(1)|temp2(2)~=1|det(a)==0

    error('矩阵或向量的大小不对应或矩阵为奇异矩阵')

    return;

end

%%化为三角矩阵过程

leng=temp1(1);N=-tril(a,-1)-triu(a,1);M=a+N;

t=1;

for i=1:leng

    t=t*M(i,i);

end

if t==0

    error('对角元素都为零！')

end

x=zeros(leng,1);

t=0;x1=zeros(leng,1);

for i=1:leng

    for j=1:leng

        if j

            x(i)=N(i,j)*x(j)+x(i);

        end

        if j>i

            x(i)=N(i,j)*x1(j)+x(i);

        end

    end

    x(i)=(x(i)+b(i))/M(i,i);

end

for i=1:leng  

    t=t+(x1(i)-x(i))^2;

end

e=E^2;

n=1;

while t>e

   t=0;

   x1=x;

   x=zeros(leng,1);

   for i=1:leng

       for j=1:leng

           if j

              x(i)=N(i,j)*x(j)+x(i);

           end

           if j>i

              x(i)=N(i,j)*x1(j)+x(i);

           end

       end

       x(i)=(x(i)+b(i))/M(i,i);

       t=t+(x1(i)-x(i))^2;

   end

   n=n+1;

   if t>1000

       break;

       error('不收敛')

   end

end

控制命令窗口命令符为：

a=[8 -3 2;4 11 -1; 6 3 12];b=[20 33 36]';

>> [x,n]=gauss_seidel(a,b,0.001)

x =
    2.9998

    2.0001

    1.0001

n =
     5

用SOR迭代法解法方程①

SOR.m文件为

function [x,n]=SOR(a,b,E,p)

temp1=size(a);temp2=size(b);

if temp1(1)~=temp1(2)|temp1(1)~=temp2(1)|temp2(2)~=1|det(a)==0

    error('矩阵或向量的大小不对应或矩阵为奇异矩阵')

    return;

end

%%化为三角矩阵过程

leng=temp1(1);N=-tril(a,-1)-triu(a,1);M=a+N;N=-a;

t=1;

for i=1:leng

    t=t*M(i,i);

end

if t==0

    error('对角元素都为零！')

end

x=zeros(leng,1);

t=0;x1=zeros(leng,1);

for i=1:leng

       for j=1:leng

        if j

            x(i)=N(i,j)*x(j)+x(i);

        end

        if j>=i

           x(i)=N(i,j)*x1(j)+x(i);

        end

    end

    x(i)=x1(i)+(x(i)+b(i))/M(i,i)*p;

end

for i=1:leng  

    t=t+(x1(i)-x(i))^2;

end

e=E^2;

n=1;

while t>e

   t=0;

   x1=x;

   x=zeros(leng,1);

   for i=1:leng

       temp=0;

       for j=1:leng

           if j

              x(i)=N(i,j)*x(j)+x(i);

           end

           if j>=i

              x(i)=N(i,j)*x1(j)+x(i);

           end

       end

       x(i)=x1(i)+(x(i)+b(i))/M(i,i)*p;

       t=t+(x1(i)-x(i))^2;

   end

   n=n+1;

   if t>1000

       break;

       error('不收敛')

   end

end

控制命令窗口命令为：

当松弛因子p=1.3时

[x,n]=SOR(a,b,0.001,1.3)

x =
    3.0000

    2.0000

    1.0001
n =

    11

当松弛因子p=1时也即是高斯-塞德尔方法

>> [x,n]=SOR(a,b,0.001,1)
x =
    2.9998

    2.0001

    1.0001

n =
     5

当松弛因子p=0.9时

[x,n]=SOR(a,b,0.001,0.9)
x =
    3.0001

    1.9999

    0.9999
n =
     6

所以当p=1时松弛因子比较好

通过比较发现 高斯-塞德尔迭代法收敛的速度明显比雅各比收敛速度快

3 用直接法解地界稠密矩阵是占优势，用迭代法解大型稀疏矩阵占优势