求解线性方程组的方法Matlab程序

求解线性方程组的方法

目录

一，线性方程组的直接解法

1.高斯消元法
(1).原理和算法
(2).Matlab函数
(3).运行实例

2.用主元法进行高斯消元
(1).原理和算法
(2).Matlab函数
(3).运行实例

3.用按比例主元法进行高斯消元
(1).原理和算法
(2).Matlab函数
(3).运行实例

二，用迭代法求解线性方程组

1.雅可比迭代法
(1).原理和算法
(2).Matlab函数
(3).运行实例

2.高斯塞达尔迭代法
(1).原理和算法
(2).Matlab函数
(3).运行实例

一，线性方程组的直接解法
1.高斯消元法
(1).原理和算法

(2).Matlab函数

function x=Gauss(A,b)
n=size(A,1);
x=zeros(n,1);
for i=1:n-1 
    for t=i+1:n%用第i行的从第i+1行开始消去，一直消到最后一行
       m(t,i)=-A(t,i)/A(i,i); %后一行首位/前一行首位得到消元乘数
       for j=i:n %第i行从第i个开始乘乘数，消去第i+1行第i个元素
           A(t,j)=A(t,j)+m(t,i)*A(i,j);
       end
       b(t)=b(t)+m(t,i)*b(i); %该行b也随之变化
    end
end%完成所有行消元，因为最后一行最后一列元素不用消，所以循环到n-1
for i=1:n
    y(i)=b(i);
end 
for i=n:-1:1
    t=n;
    while i<t
          y(i)=y(i)-A(i,t)*x(t);
          t=t-1;
    end
    x(i)=y(i)/A(i,i);
end%回代求解

(3).运行实例
求解线性方程组：
x1+ x2 + x3+x4=4
2x1+ x2 + x3+x4=5
3x1+2x2+x3 +x4=7
4x1+3x2+2x3+x4=10

2.用主元法进行高斯消元
(1).原理和算法

(2).Matlab函数

function x=Gauss(A,b)
n=size(A,1);
x=zeros(n,1);
for i=1:n-1 
   q=i;%q用来存储列元素最大的列数，假设第i行最大
    for p=i:n%第i列内从第i行开始寻找最大元素
        if A(q,i)<A(p,i)
           q=p; 
        end
    end
    A([i q],:)=A([q i],:); %元素最大行与第i行交换
    temp=b(i);b(i)=b(q);b(q)=temp; %该两行b也随之交换
for t=i+1:n %用第i行的从第i+1行开始消去，一直消到最后一行
       m(t,i)=-A(t,i)/A(i,i); %后一行首位/前一行首位得到消元乘数
       for j=i:n %第i行从第i个开始乘乘数，消去第i+1行第i个元素
           A(t,j)=A(t,j)+m(t,i)*A(i,j);
       end
       b(t)=b(t)+m(t,i)*b(i); %该行b也随之变化
    end
end%完成所有行消元，因为最后一行最后一列元素不用消，所以循环到n-1
for i=1:n
    y(i)=b(i);
end 
for i=n:-1:1
    t=n;
    while i<t
          y(i)=y(i)-A(i,t)*x(t);
          t=t-1;
    end
    x(i)=y(i)/A(i,i);
end%回代求解

(3).运行实例
求解线性方程组：
x1+ x2 + x3+x4=4
2x1+ x2 + x3+x4=5
3x1+2x2+x3 +x4=7
4x1+3x2+2x3+x4=10

3.用按比例主元法进行高斯消元
(1).原理和算法

(2).Matlab函数

function x=Gauss(A,b)
n=size(A,1);
x=zeros(n,1);
for i=1:n-1 

for r=i:n %选出每行的最大元素并和该行第一个元素相除
        h=r; %h用来存储行元素最大的行数，假设第i个最大
        for k=r:n%第r行内从第1个开始寻找最大元素
           if A(r,h)<A(r,k)
             h=k; 
           end%最终h保留的为元素最大的列数
        end
        D(r)=A(r,r)/A(r,h);%存储在数列内
    end
    q=i;%q用来存储列元素最大的列数，假设第i行最大
    for p=i:n%第i列内从第i行开始寻找最大比值
       if D(q)<A(p)
         q=p; 
       end%最终q保留的为元素最大的行数
    end
    A([i q],:)=A([q i],:);%元素最大行与第i行交换
    temp=b(i);b(i)=b(q);b(q)=temp; %该两行b也随之交换
for t=i+1:n%用第i行的从第i+1行开始消去，一直消到最后一行
       m(t,i)=-A(t,i)/A(i,i); %后一行首位/前一行首位得到消元乘数
       for j=i:n %第i行从第i个开始乘乘数，消去第i+1行第i个元素
           A(t,j)=A(t,j)+m(t,i)*A(i,j);
       end
       b(t)=b(t)+m(t,i)*b(i); %该行b也随之变化
    end
end%完成所有行消元，因为最后一行最后一列元素不用消，所以循环到n-1
for i=1:n
    y(i)=b(i);
end 
for i=n:-1:1
    t=n;
    while i<t
          y(i)=y(i)-A(i,t)*x(t);
          t=t-1;
    end
    x(i)=y(i)/A(i,i);
end%回代求解

(3).运行实例
求解线性方程组：
x1+ x2 + x3+x4=4
2x1+ x2 + x3+x4=5
3x1+2x2+x3 +x4=7
4x1+3x2+2x3+x4=10

二，用迭代法求解线性方程组
1.雅可比迭代法
(1).原理和算法

(2).Matlab函数

function x=jacobi(A,b,x0,tol,max) 
[n,m]=size(A);
xold=x0;
C=-A;
for i=1:n 
    C(i,i)=0;
end 
for i=1:n 
    C(i,:)=C(i,:)/A(i,i); 
end
for i=1:n 
    d(i,1)=b(i)/A(i,i); 
end
i=1;
while i<=max 
      xnew=C*xold+d; %雅可比迭代直接使用上一次的解
      E=norm(xnew-xold)
      if norm(xnew-xold)<=tol 
         x=xnew;
         disp('Accuracy achieved!') 
         return; 
      else 
         xold=xnew; 
      end 
      disp([xnew']); 
      i=i+1; 
end 
x=xnew; 

(3).运行实例

运行：

>> A=[3,-1,1;3,6,2;3,3,7];
b=[1,0,4]';
tol=0.0001;
max=20;
x0=[0;0;0;];
>> x=jacobi(A,b,x0,tol,max)
E =

    0.6615    0.3333         0    0.5714


E =

    0.4292    0.1429   -0.3571    0.4286


E =

    0.2839    0.0714   -0.2143    0.6633


E =

    0.0607    0.0408   -0.2568    0.6327


E =

    0.0406    0.0368   -0.2313    0.6640


E =

    0.0127    0.0349   -0.2398    0.6548


E =

    0.0060    0.0352   -0.2357    0.6592


E =

    0.0025    0.0350   -0.2373    0.6574


E =

    0.0010    0.0351   -0.2366    0.6581


E =

   4.3591e-04    0.0351   -0.2369    0.6578


E =

   1.7968e-04    0.0351   -0.2368    0.6579


E =7.6131e-05

Accuracy achieved!
n=12
x =

    0.0351
   -0.2369
    0.6579

2.高斯塞达尔迭代法
(1).原理和算法

(2).Matlab函数

function x=gs(A,b,x0,tol,max) 
[n,m]=size(A);
xold=x0;
C=-A;
for i=1:n 
    C(i,i)=0;
end 
for i=1:n 
    C(i,:)=C(i,:)/A(i,i); 
end
for i=1:n 
    d(i,1)=b(i)/A(i,i); 
end
i=1;
while i<=max 
      xnew=xold;
      for j=1:n
          xnew(j,1)=C(j,:)*xnew+d(j,1); %高斯塞达尔迭代保持解的实时更新
      end
      E=norm(xnew-xold)
      if norm(xnew-xold)<=tol 
         x=xnew;
         disp('Accuracy achieved!') 
         return; 
      else 
         xold=xnew; 
      end 
      disp([xnew']); 
      i=i+1; 
end
x=xnew;

(3).运行实例

运行：

>> A=[3,-1,1;3,6,2;3,3,7];
b=[1,0,4]';
tol=0.0001;
max=20;
x0=[0;0;0;];
>> x=gs(A,b,x0,tol,max)

E =

0.6236 0.3333   -0.1667    0.5000

E =
0.2582   0.1111   -0.2222    0.6190

E =
0.0661  0.0529   -0.2328    0.6485

E =
0.0154  0.0396   -0.2360    0.6556


E =

    0.0039  0.0361   -0.2366    0.6573


E =

   9.1975e-04

    0.0354   -0.2368    0.6578


E =

   2.2862e-04

    0.0352   -0.2368    0.6579


E =

   5.4682e-05

Accuracy achieved!
n=8
x =

    0.0351
   -0.2368
    0.6579