第3 章方程组的数值解法

zaiziwamei

于 2022-10-06 11:10:43 发布

阅读量1.6k

点赞数 2

分类专栏：数值计算方法

本文链接：https://blog.csdn.net/qq_40940944/article/details/127175796

版权

数值计算方法专栏收录该内容

1 篇文章 1 订阅

订阅专栏

系列文章目录

前言

本文主要介绍方程组的一些数值解法。

3.1 Gauss消去法

3.1.1 基本思想

Gauss消去法的基本思想是用原方程组中的一些方程线性组合出一个新的三角形方程组，然后回代参数求解方程。

3.1.2 求解过程

参见褚衍东《数值计算方法》P49-51。

3.1. 3 实例求解

用MATLAB语言求解如下方程组：
在这里插入图片描述在MATLAB程序中新建Gauss.m文件并编写以下代码

function f = Gauss(A,b)
    % 函数功能：顺序高斯消去法解线性方程组Ax=0
    % f = Gauss(A,b),A为系数矩阵，b为右端列向量，x为解向量。
    A=[A';b]';%先将矩阵A求转置，之后将向量b连接到转置矩阵A的下面，最后再对连接后的矩阵呢个求转置。
    n=length(b);%获取b的长度
    % 1消去过程
    for k =1:n-1
        for i=k+1:n
            m =A(i,k)/A(k,k);%m为课本中l
            %fprintf('m%d%d=%f\n',i,k,m);
            for j=k:n+1
                A(i,j)=A(i,j)-m*A(k,j);%
            end
        end
        %fprintf('A%d=',k+1);% 输出每对一列进行置0变换后的矩阵
        %disp(A) % 输出矩阵A（每对一列进行置0变换后的矩阵）
    end
    %A
    % 2回带过程
        % 在求解过程中将Xi,将Xi放在bi的位置。
        % 回带过程中系数矩阵不在改变。
    A(n,n+1)=A(n,n+1)/A(n,n);%求解Xn
    for i=n-1:-1:1
        s=0;
        for j=i+1:n
            s=s+A(i,j)*A(j,n+1);
        end
        A(i,n+1)=(A(i,n+1)-s)/A(i,i);%求解Xi
    end
    A(:,n+1)%输出方程的解
end

在MATLAB命令窗口输入：

A=[1 1 0 -4;-1 1 1 3;1 3 5 -4;0 1 2 -1];
b=[1 -2 -4 -2];
Gauss(A,b)

输出：

ans =

     1
     0
    -1
     0

即方程的解:x₁=1,x₂=0;x₃=-1;x₄=0.

3.2 列主元Gauss消去法

3.2.1 基本思想

在进行第 k(k=1,2,…,n-1)步消元时，从第k列的 a_kk及其以下的各元素中选取绝对值最大的元素，然后通过行变换将它交换到主元素a_kk的位置上，再进行消元。

3.2.2 实例求解

用MATLAB求解以下方程组：
在这里插入图片描述

在MATLAB程序中新建Gauss_lzy.m文件并编写以下代码

function [x,det,flag] = Gauss_lzy(A,b)
% 函数功能：求线性方程组的列主元Gauss消去法
% A为线性方程组的系数矩阵，b为线性方程组的右端项，x为方程组的解，det为系数矩阵A的行列式
% flag为指标向量，flag='failure'表示计算失败，flag='OK'表示计算成功

[n,m]=size(A);nb=length(b);
% 当方程组行与列的维数不相等时，停止计算并输出出错信息
if n~=m
    error('The rows and columns of matrix A must be equal!');
    return;
end
% 当方程组与右端的维数不匹配时，停止计算，并输出错误信息
if m~=nb
    error('The colums if A must be equal the length of b!');
    return;
end
% 开始计算，先赋初值
flag ='OK';det=1;x=zeros(n,1);
for k=1:n-1
    % 选主元
    max1=0;
    for i =k:n
        if abs(A(i,k))>max1
            max1=abs(A(i,k));r=1;
        end
    end
    if max1<1e-10
        flag='failure';return;
    end
    %交换两行
    if r>k
        for j=k:n
            z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;
        end
        z=b(k);b(k)=b(r);b(r)=z;det=-det;
    end
    % 消元过程
    for i=k+1:n
        m=A(i,k)/A(k,k);
        for j=k+1:n
            A(i,j)=A(i,j)-m*A(k,j);
        end
        b(i)=b(i)-m*b(k);
    end
    det=det*A(k,k);
end
det=det*A(n,n);
%回带过程
if abs(A(n,n))<1e-10
    flag='failure';
    return;
end
for k=n:-1:1
    for j=k+1:n
        b(k)=b(k)-A(k,j)*x(j);
    end
    x(k)=b(k)/A(k,k);
end
x(k)=b(k)/A(k,k);
end

在MATLAB命令窗口输入：

A=[0.000100 1;1 1];
b=[1 2];
[x,det,flag] = Gauss_lzy(A,b)

输出：

x =

    1.0001
    0.9999


det =

   -0.9999


flag =

    'OK'

即方程的解:x₁=1,x₂=0.9999.

3.3 矩阵的三角分解法

3.3.1 基本思想

设方程组的矩阵形式为Ax=b，三角分解法就是将系数矩阵A分解为一个下三角矩阵L和一个上三角矩阵U之积：A=LU，然后依次解两个三角形方程组Ly=b和Ux=y，而得到原方程组的解。

3.3.2 实例求解

用MATLAB求解以下方程组：
在这里插入图片描述

在MATLAB程序中新建LU_x.m文件并编写以下代码

function [y,x] = LU_x(A,b)
% 函数功能：选主元LU分解法解方程
% f = f[y,x] = LU_x(A,b),A为系数矩阵，b为右端列向量，x为解向量。
b=b';A=[A',b]';n=length(b');x=zeros(n,1);y=zeros(n,1);
U=zeros(n);L=zeros(n);
for k=1:n
    U(1,k)=A(1,k);
    L(k,1)=A(k,1)/U(1,1);
end

for i=2:n
    for j=i:n
        lu=0;
        lu1=0;
        for k=1:i-1
            lu=lu+L(i,k)*U(k,j);
            lu1=lu1+L(i,k)*U(k,j);
        end
        U(i,j)=A(i,j)-lu;
        L(j,i)=(A(j,i)-lu1)/U(i,i);
    end 
end
%L
%U
%求解方程Ly=b
ly=0;
% y(1)=b(1);
for i=1:n
    ly=0;
    for j=1:i-1
        ly =ly+L(i,j)*y(j);
    end
    y(i)=b(i)-ly;
end
%求解方程Ux=y
for i=n:-1:1
    ly1=0;
    for j=i+1:n
        ly1=ly1+U(i,j)*x(j);
    end
    x(i)=(y(i)-ly1)/U(i,i);
end
end

在MATLAB命令窗口输入：

A=[2 4 2 6;4 9 6 15;2 6 9 18; 6 15 18 40];
b=[9 23 22 47];
[y,x]=solve_equction.LU_x(A,b)

即方程的解:x₁=0.5,x₂=2，x₃=3，x₄=-1.

3.4 追赶法

3.4.1 基本思想

方程组的矩阵形式为Ax=b，本文中A为三对角矩阵，追赶法就是将系数矩阵A分解为一个下二对角矩阵L和一个上二对角矩阵U之积：A=LU，然后依次解两个三角形方程组Ly=b和Ux=y，而得到原方程组的解。

3.4.2 实例求解

用MATLAB求解以下方程组：
在这里插入图片描述

在MATLAB程序中新建chase.m文件并编写以下代码

function f=chase(a,b,c,d)
%clear all;
% a=[2,2,2,2];% 将三对角元素存入向量a,b,c
% b=[-1,-1,-1];
% c=[-1,-1,-1];
% d=[1,0,0,1];% 将右端向量存入r
n=length(a);
b=[0,b];% 由于b的序号从2开始，故前面需要多增加一项
% 用追赶法计算u(1).v(1),……u(n-1),v(n-1),u(n)
u(1)=d(1)/a(1);
v(1)=c(1)/a(1);
for k=2:n-1
    u(k)=(d(k)-u(k-1)*b(k))/(a(k)-v(k-1)*b(k));
    v(k)=c(k)/(a(k)-v(k-1)*b(k));
end
u(n)=(d(n)-u(n-1)*b(n))/(a(n)-v(n-1)*b(n));
%用追赶法计算x(n),x(n-1),……，x(1)
x(n)=u(n);
for k=n-1:-1:1
    x(k)=u(k)-v(k)*x(k+1);
end
fprintf('三队方程组的解为\n')
for k=1:n
    fprintf('x(%1d)=%10.8f\n',k,x(k));
end
end

输出：

三队方程组的解为
x(1)=1.00000000
x(2)=1.00000000
x(3)=1.00000000
x(4)=1.00000000

即方程的解:x₁=1,x₂=1，x₃=1，x₄=1.

3.5 平方根法

3.5.1 基本思想

平方根分解矩阵是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积。

3.5.2 实例求解

用MATLAB求解以下方程组：
在这里插入图片描述本题只涉及如何将求其下三角矩阵L，至于后续方程的求解与前文矩阵的三角分解法类似，这里不再求解。
在MATLAB程序中新建Cholesky.m文件并编写以下代码

function L=Cholesky(A)
% 功能：(平方根法)正定对称矩阵LU分解的Cholesky法
% A为LU分解矩阵，L为下三角矩阵
n=length(A);L=zeros(n);

for j=1:n
    delta=A(j,j);
    for k=1:j-1
        delta=delta-(L(j,k)^2);
    end
    L(j,j)=sqrt(delta);
    for i=j+1:n
        L(i,j)=A(i,j);
        for k=1:j-1
            L(i,j)=L(i,j)-L(i,k)*L(j,k);
        end
        L(i,j)=L(i,j)/L(j,j);
    end
end
end

在MATLAB命令框中输入

A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10]
Cholesky(A)

输出：
ans =

    2.2361         0         0         0
    3.1305    0.4472         0         0
    2.6833   -0.8944    1.4142         0
    2.2361         0    2.1213    0.7071

3.6 改进的平方根法

3.6.1 基本思想

把一个对称正定的矩阵A做LDL^T分解使得L中主对角线元素均为1，D中除主对角线以外的元素均为0.

3.6.2 实例求解

用MATLAB求解以下方程组：
在这里插入图片描述本题只涉及如何将求其下三角矩阵L，对角矩阵D，至于后续方程的求解与前文矩阵的三角分解法类似，这里不再求解。
在MATLAB程序中新建LDL.m文件并编写以下代码

function [L,D]=LDL(A)
% 功能：（改进的平方根法）正定对称矩阵LU分解的Cholesky法，
% A为要LU分解的矩阵，L为下三角矩阵，D为对角矩阵
n=length(A);L=zeros(n);d=zeros(1,n);T=zeros(n);
for k=1:n
    d(k)=A(k,k);
    for j=1:k-1
        d(k)=d(k)-L(k,j)*T(k,j);
    end
    for i=k+1:n
        T(i,k)=A(i,k);
        for j=1:k-1
            T(i,k)=T(i,k)-T(i,j)*L(k,j);
        end
        L(i,k)=T(i,k)/d(k);
    end
    L(k,k)=1;
end
D=diag(d);
end

在MATLAB命令框中输入

A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10];
[L,D]=LDL(A)

输出：
L =

    1.0000         0         0         0
    1.4000    1.0000         0         0`在这里插入代码片`
    1.2000   -2.0000    1.0000         0
    1.0000         0    1.5000    1.0000


D =

    5.0000         0         0         0
         0    0.2000         0         0
         0         0    2.0000         0
         0         0         0    0.5000

3.7 函数整理与调用

将这些函数整理在solve_fun.m文件中，在solve_main.m文件中根据需要调用。
solve_fun.m中代码如下：

function fun_equation=solve_fun
    fun_equation.Gauss = @Gauss;
    fun_equation.Gauss_lzy = @Gauss_lzy;
    fun_equation.GaussXQAllMain = @GaussXQAllMain;
    fun_equation.LU_x = @LU_x;
    fun_equation.chase = @chase;
    fun_equation.Cholesky = @Cholesky;
    fun_equation.LDL = @LDL;
end

function f = Gauss(A,b)
    % 函数功能：顺序高斯消去法解线性方程组Ax=0
    % f = Gauss(A,b),A为系数矩阵，b为右端列向量，x为解向量。
    A=[A';b]';%先将矩阵A求转置，之后将向量b连接到转置矩阵A的下面，最后再对连接后的矩阵呢个求转置。
    n=length(b);%获取b的长度
    % 1消去过程
    for k =1:n-1
        for i=k+1:n
            m =A(i,k)/A(k,k);%m为课本中l
            %fprintf('m%d%d=%f\n',i,k,m);
            for j=k:n+1
                A(i,j)=A(i,j)-m*A(k,j);%
            end
        end
        %fprintf('A%d=',k+1);% 输出每对一列进行置0变换后的矩阵
        %disp(A) % 输出矩阵A（每对一列进行置0变换后的矩阵）
    end
    %A
    % 2回带过程
        % 在求解过程中将Xi,将Xi放在bi的位置。
        % 回带过程中系数矩阵不在改变。
    A(n,n+1)=A(n,n+1)/A(n,n);%求解Xn
    for i=n-1:-1:1
        s=0;
        for j=i+1:n
            s=s+A(i,j)*A(j,n+1);
        end
        A(i,n+1)=(A(i,n+1)-s)/A(i,i);%求解Xi
    end
    A(:,n+1)%输出方程的解
end
function [x,det,flag] = Gauss_lzy(A,b)
% 函数功能：求线性方程组的列主元Gauss消去法
% A为线性方程组的系数矩阵，b为线性方程组的右端项，x为方程组的解，det为系数矩阵A的行列式
% flag为指标向量，flag='failure'表示计算失败，flag='OK'表示计算成功

[n,m]=size(A);nb=length(b);
% 当方程组行与列的维数不相等时，停止计算并输出出错信息
if n~=m
    error('The rows and columns of matrix A must be equal!');
    return;
end
% 当方程组与右端的维数不匹配时，停止计算，并输出错误信息
if m~=nb
    error('The colums if A must be equal the length of b!');
    return;
end
% 开始计算，先赋初值
flag ='OK';det=1;x=zeros(n,1);
for k=1:n-1
    % 选主元
    max1=0;
    for i =k:n
        if abs(A(i,k))>max1
            max1=abs(A(i,k));r=1;
        end
    end
    if max1<1e-10
        flag='failure';return;
    end
    %交换两行
    if r>k
        for j=k:n
            z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;
        end
        z=b(k);b(k)=b(r);b(r)=z;det=-det;
    end
    % 消元过程
    for i=k+1:n
        m=A(i,k)/A(k,k);
        for j=k+1:n
            A(i,j)=A(i,j)-m*A(k,j);
        end
        b(i)=b(i)-m*b(k);
    end
    det=det*A(k,k);
end
det=det*A(n,n);
%回带过程
if abs(A(n,n))<1e-10
    flag='failure';
    return;
end
for k=n:-1:1
    for j=k+1:n
        b(k)=b(k)-A(k,j)*x(j);
    end
    x(k)=b(k)/A(k,k);
end
x(k)=b(k)/A(k,k);
end
function [x,XA]=GaussXQAllMain(A,b)
% 函数功能：求线性方程组的全主元Gauss消去法
% A为线性方程组的系数矩阵，b为线性方程组的常数向量，XA为消元后的系数矩阵
N=size(A);n=N(1);
idex_1=0;
index_r=0;
order=1:n;%记录未知数顺序的向量
for i=1:(n-1)
    me=max(max(abs(A(i:n,i:n))));% 选取全主元
    for k=i:n
        for r=i:n
            if(abs(A(k,r))==me)
                index_1=k;
                index_r=r;%保存主元所在的列和行
                k=n;
                break;
            end
        end
    end
end
temp=A(i,1:n);
A(i,1:n)=A(index_1,1:n);
A(index_1,1:n)=temp;
bb=b(index_1);
b(index_1)=b(i);
b(i)=bb;% 交换主行
temp=A(1:n,i);
A(1:n,i)=A(1:n,index_r);
A(1:n,index_r)=temp;% 交换主列
pos=order(i);
order(i)=order(index_r);
order(index_r)=pos;% 主列的交换会造成未知数顺序的变化
for j=(i+1):n
    if A(i,j)==0
        dis('对角线元素为0');
        return
    end
    l=A(j,i);
    m=A(i,i);
    A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m;
    b(j)=b(j)-1*b(i)/m;
end
A
b
x=b%SolveUpTriangle(A,b);
y=zeros(n,1);
for i=1:n
    for j=1:n
        if order(j)==i
            y(i)=x(j);
        end
    end
end% 恢复未知数原来的顺序
x=y
XA=A
end
function [y,x] = LU_x(A,b)
% 函数功能：选主元LU分解法解方程
% f = f[y,x] = LU_x(A,b),A为系数矩阵，b为右端列向量，x为解向量。
b=b';A=[A',b]';n=length(b');x=zeros(n,1);y=zeros(n,1);
U=zeros(n);L=zeros(n);
for k=1:n
    U(1,k)=A(1,k);
    L(k,1)=A(k,1)/U(1,1);
end

for i=2:n
    for j=i:n
        lu=0;
        lu1=0;
        for k=1:i-1
            lu=lu+L(i,k)*U(k,j);
            lu1=lu1+L(i,k)*U(k,j);
        end
        U(i,j)=A(i,j)-lu;
        L(j,i)=(A(j,i)-lu1)/U(i,i);
    end 
end
%L
%U
%求解方程Ly=b
ly=0;
% y(1)=b(1);
for i=1:n
    ly=0;
    for j=1:i-1
        ly =ly+L(i,j)*y(j);
    end
    y(i)=b(i)-ly;
end
%求解方程Ux=y
for i=n:-1:1
    ly1=0;
    for j=i+1:n
        ly1=ly1+U(i,j)*x(j);
    end
    x(i)=(y(i)-ly1)/U(i,i);
end
end
function f=chase(a,b,c,d)
%clear all;
% a=[2,2,2,2];% 将三对角元素存入向量a,b,c
% b=[-1,-1,-1];
% c=[-1,-1,-1];
% d=[1,0,0,1];% 将右端向量存入r
n=length(a);
b=[0,b];% 由于b的序号从2开始，故前面需要多增加一项
% 用追赶法计算u(1).v(1),……u(n-1),v(n-1),u(n)
u(1)=d(1)/a(1);
v(1)=c(1)/a(1);
for k=2:n-1
    u(k)=(d(k)-u(k-1)*b(k))/(a(k)-v(k-1)*b(k));
    v(k)=c(k)/(a(k)-v(k-1)*b(k));
end
u(n)=(d(n)-u(n-1)*b(n))/(a(n)-v(n-1)*b(n));
%用追赶法计算x(n),x(n-1),……，x(1)
x(n)=u(n);
for k=n-1:-1:1
    x(k)=u(k)-v(k)*x(k+1);
end
fprintf('三队方程组的解为\n')
for k=1:n
    fprintf('x(%1d)=%10.8f\n',k,x(k));
end
end
function L=Cholesky(A)
% 功能：(平方根法)正定对称矩阵LU分解的Cholesky法
% A为LU分解矩阵，L为下三角矩阵
n=length(A);L=zeros(n);

for j=1:n
    delta=A(j,j);
    for k=1:j-1
        delta=delta-(L(j,k)^2);
    end
    L(j,j)=sqrt(delta);
    for i=j+1:n
        L(i,j)=A(i,j);
        for k=1:j-1
            L(i,j)=L(i,j)-L(i,k)*L(j,k);
        end
        L(i,j)=L(i,j)/L(j,j);
    end
end
end
function [L,D]=LDL(A)
% 功能：（改进的平方根法）正定对称矩阵LU分解的Cholesky法，
% A为要LU分解的矩阵，L为下三角矩阵，D为对角矩阵
n=length(A);L=zeros(n);d=zeros(1,n);T=zeros(n);
for k=1:n
    d(k)=A(k,k);
    for j=1:k-1
        d(k)=d(k)-L(k,j)*T(k,j);
    end
    for i=k+1:n
        T(i,k)=A(i,k);
        for j=1:k-1
            T(i,k)=T(i,k)-T(i,j)*L(k,j);
        end
        L(i,k)=T(i,k)/d(k);
    end
    L(k,k)=1;
end
D=diag(d);
end

solve_main.m文件中代码如下：
在调用某个方法时仅需要去掉其前面的注释。

solve_equction= solve_fun;
%% Gauss消去法
%A=[1 1 0 -4;-1 1 1 3;1 3 5 -4;0 1 2 -1];
%b=[1 -2 -4 -2];
%solve_equction.Gauss(A,b)
%% 列主元Gauss消去法
%A=[0.000100 1;1 1];
%b=[1 2];
%[x,det,flag] = solve_equction.Gauss_lzy(A,b) 
%% 全主元Gauss消去法
%A=[0.000100 1;1 1];
%b=[1 2];
%[x,XA]=solve_equction.GaussXQAllMain(A,b)
%% 选主元LU分解法
%A=[2 4 2 6;4 9 6 15;2 6 9 18; 6 15 18 40];
%b=[9 23 22 47];
%[y,x]=solve_equction.LU_x(A,b)
%% 追赶法求解
% a=[2,2,2,2];% 将三对角元素存入向量a,b,c
% b=[-1,-1,-1];
% c=[-1,-1,-1];
% d=[1,0,0,1];% 将右端向量存入r
% solve_equction.chase(a,b,c,d)
%% 平方根法求解
%A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10];
%solve_equction.Cholesky(A)
%% 改进的平方根法求解
%A=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10];
%[L,D]=solve_equction.LDL(A)