数值分析中求解方程组或方程组的根的基本算法（Matlab实现）

0 方程

A1 = [5 6 7
     2 2 3
     3 5 1];
b1 = [12 9 15];
% x1=33 x2=-15 x3=-9 
% 三元一次方程组在线计算器：https://www.osgeo.cn/app/sc290

A2 = [2 3 4 -5
    6 7 -8 9
    10 11 12 13
    14 15 16 17];
b2 = [-6 96 312 416];
% x1=3 x2=5 x3=7 x4=11 

% 四元一次方程组在线计算器：https://www.osgeo.cn/app/sc291

1. 利用高斯消元法求解线性方程组。

            x1  +  x2 +  x3  + x4 = 4

           2x1  +  x2 +  x3  + x4 = 5

           3x1  + 2x2 +  x3  + x4 = 7

           4x1  + 3x2 + 2x3 + x4 = 10

运行结果：

代码：

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%回代求解

---

2. 利用雅可比迭代法求解线性方程组。

           3x1 - x2  +  x3 = 1
           
           3x1 + 6x2 + 2x3 = 0
           
           3x1 + 3x2 + 7x3 = 4

运行结果：

代码：

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. 利用牛顿迭代法求方程的根。

          x * e^x - 1 = 0

运行结果：

代码：

function [X_k,x0,counter]=newton(a,err,f_x)
counter=0;%计数器初始化
X_k=(0);%序列初始化
while(sign(f_x(a))==0)%如果a的函数值为0停止迭代输出x0为a的值，迭代次数为0
    x0=a;
    X_k(1)=a;
    counter=0;
    return;
end
while(sign(f_x(a))~=0&&abs(f_x(a)-a)>=err)%a的函数不为零且不满足误差精度进行迭代
    X_k(counter+1)=a;%将a赋给序列第一个值
    a=f_x(a);%迭代
    counter=counter+1;%每迭代一次计数器加一
    X_k(counter+1)=a;%将当前值赋值给序列
    x0=a;%最终输出最后一个迭代值作为根
    newton(a,err,f_x);%不满足精度时重复调用迭代函数本身进行迭代
end
End

---

4. 利用二分法求方程的根。

          y = x^2 + x - 1

运行结果：

代码：

function root=bisect(fun,a,b,eps)
n=1+round((log(b-a)-log(eps))/log(2));
fa=feval(fun,a);fb=feval(fun,b);
for i=1:n
    c=(b+a)/2;
    fc=feval(fun,c);
    if fc*fa<0
        b=c;fb=fc;
    else
        a=c;fa=fc;
    end
end
root=c;

---

5. 利用LU分解法求解线性方程组。

          3x1 - x2  +  x3 = 1
          
          3x1 + 6x2 + 2x3 = 0
          
          3x1 + 3x2 + 7x3 = 4

运行结果：

代码：

function x=LUfenjiefa(A,b)
n=length(b);
k=2;
X=A;
Y=b;
U(1,1:n)=A(1,1:n);
L(2:n,1)=A(2:n,1)/U(1,1);
for k=2:n
    U(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);
    L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k))/U(k,k);
end
L
U
%用向前消去法解下三角方程组Ly=b
y=zeros(n,1);
y(1)=b(1);
for k=2:n
    y(k)=b(k)-L(k,1:k-1)*y(1:k-1);
end
y
%用回代法解上上角方程组Ux=y
x=zeros(n,1);
x(n)=y(n)/U(n,n);
for k=n-1:-1:1
    x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);
end
end

6 选取列主元的Doolittle分解法求解线性方程组

function x = Doolittle(A,b)
    % 假设A的各阶顺序主子式均不为零
    n = length(A);
    L = eye(n);
    U = zeros(n);
    %length(A)返回矩阵A最大维数；
    %eye(n) 返回一个主对角线元素为 1 且其他位置元素为 0 的 n×n 单位矩阵。
    
    for k = 1:n-1 
        for i = k:n
            % 计算下三角矩阵L第k列元素的一部分
            s(i) = A(i,k)-L(i,1:k-1)*U(1:k-1,k);
        end

        %选列主元
        [~,q] = max(abs(s(k:n)));
        q = q+k-1; %还原s(i)的索引下标
        % 若q!=k
        if q>k %绝对值最大的列主元不在第q行即当前第一行
        % 交换A的第k行与第q行
            t_A = A(k,:);
            A(k,:) = A(q,:);
            A(q,:) = t_A;
        % 交换下三角矩阵L第k列元素s_k与s_q（l_kk与l_qk）
            t_s = s(k);
            s(k) = s(q);
            s(q) = t_s;
        % 将下三角矩阵L的第k行前(k-1)个元素与第q行的前(k-1)个元素交换
            t_L = L(k,1:k-1);
            L(k,1:k-1) = L(q,1:k-1);
            L(q,1:k-1) = t_L;
        % 交换b_k与b_q
            t_b = b(k);
            b(k) = b(q);
            b(q) = t_b;
        end
        % 计算矩阵U的第k行元素及矩阵L的第k列元素
        U(k,k) = s(k);
        for i = k+1:n
            U(k,i) = A(k,i)-L(k,1:k-1)*U(1:k-1,i);
            L(i,k) = s(i)/U(k,k);
        end
    end
    % 求解上三角矩阵第n行的唯一待定元素u_nn
    U(n,n) = A(n,n)-L(n,1:n-1)*U(1:n-1,n);
    % 解单位下三角方程组Ly=b
    y = zeros(n,1);
    for k = 1:n
        if k==1 
            y(1) = b(1);
        else 
            y(k) = b(k)-L(k,1:k-1)*y(1:k-1);
        end
    end
    % 解上三角方程组Ux=y
    x = zeros(n,1);
    for k = n:-1:1
        if k==n
            x(n) = y(n)/U(n,n);
        else
            x(k) = (y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);
        end
    end
    x = x';
end 

7 解线性方程组系数矩阵为三对角矩阵的追赶法

%行对角占优条件：
%   (1)比较绝对值，第一行的主对角元大于第一行的另一个非零元，且绝对值均不为零
%   (2)比较绝对值，最后一行的主对角元大于最后一行的另一个非零元，且绝对值均不为零
%   (3)第一行和最后一行之外的其余各行，主对角元的绝对值大于另外两个非零元的绝对值之和且均不为0
 
%说明：
%   矩阵A的主对角线元素为d(1)~d(n)
%   次下对角线元素为a(2)~a(n)
%   次上对角线元素为c(1)~c(n-1)
 
function x = ForwardBackward(a,b,c,d)
    n = length(d);
    % 输入4个一维数组a,c,d(存放A的三条对角线上的元素),b(Ax=b，列向量)
    % 若d(1)==0，则输出信息"方法失败！"，算法终止。
    if d(1)==0
        warning('方法失败！')
        return
    end

    p(1) = d(1);q(1)=c(1)/d(1);
    for i = 2:n-1
        p(i) = d(i)-a(i)*q(i-1);
        % 若p(i)=0,则输出信息"方法失败！"，算法终止(可以证明系数矩阵A满足对角占优条件，则p(1)~p(n)均不为0)
        if p(i)==0
            warning('方法失败！')
        return
        end
        q(i) = c(i)/p(i);
    end
    p(n) = d(n)-a(n)*q(n-1);

    % 算法终止条件判断
    if p(n)==0
        warning('方法失败！')
        return
    end

    % 解下三角方程组Ly=b
    y(1)=b(1)/p(1);
    for i = 2:n
        y(i) = (b(i)-a(i)*y(i-1))/p(i);
    end

    % 解上三角方程组Ux=y
    x(n) = y(n);
    for i = n-1:-1:1
        x(i) = y(i)-q(i)*x(i+1);
    end
end

8 利用SOR(逐次超松弛)迭代法解线性方程组

%% 用SOR(Successive Over-Relaxation,逐次超松弛)迭代法解n阶线性方程组Ax=b
function x = SOR(A,b,omega,x0,eps,N)
    % 其中松弛因子omega(ω),初始向量x0，精度eps,最大迭代次数N
     
    %步骤1：输入矩阵A，右端向量b,松弛因子ω和初始向量x0,精度eps和最大迭代次数N。令k=0
    n=length(x0);
    x=ones(n,1);k=0;
    %x其实可以设定为任意值，不影响后面的迭代和覆盖过程，设为全部元素为1的列向量纯粹为了表示方便
     
    %步骤2：当k≤N时，执行步骤2.1~2.3
    while k<=N
        %步骤2.1：对i=1,2,...n，计算x_i(k+1)表达式
        for i=1:n
            x(i) = omega*(-A(i,1:i-1)* x(1:i-1)-A(i,i+1:n)*x0(i+1:n)+b(i))/A(i,i)+(1-omega)*x0(i);
        end
        %步骤2.2：令k:=k+1
        k=k+1;
        %步骤2.3：若||x(k+1)-x(k)||<eps，则算法停止，输出方程组近似解x(k+1);否则令x_k:=x_k+1。
        if(norm(x-x0,inf)<eps),break;end
        x0=x;
    end
     
    %步骤3：输出信息“算法超出最大迭代次数！”，算法终止。
    if k>N
        warning('算法超出最大迭代次数！');
    else
        disp(['迭代次数= ',num2str(k)]);
        x
    end
end