MATLAB求解线性方程组的八种方法

2023-04-24 00:59| 来源: 网络整理| 查看: 265

MATLAB求解线性方程组的八种方法

求解线性方程分为两种方法–直接法和迭代法常见的方法一共有8种直接法 Gauss消去法 Cholesky分解法迭代法 Jacobi迭代法 Gauss-Seidel迭代法超松弛迭代法共轭梯度法 Bicg迭代法 Bicgstab迭代法

这里我就从计算代码的角度来讲解，在下面也会按照上面这个顺序给出代码，遇到方程组直接带入已知条件就可以得到答案。

适用条件 Gauss消去法：求解中小规模线性方程（阶数不过1000），一般用于求系数矩阵稠密而且没有任何特殊结构的线性方程组

Cholesky分解法：对称正定方程优先使用，系数矩阵A是n阶对称正定矩阵

Jacobi迭代法非奇异线性方程组，分量的计算顺序没有关系

Gauss-Seidel迭代法与Jacobi迭代法相似，但计算的分量不能改变

超松弛迭代法Jacobi迭代法和Gauss-Seidel迭代法的加速版，由Gauss-Seidel迭代法改进而来，速度较快

共轭梯度法需要确定松弛参数w,只有系数矩阵具有较好的性质时才可以找到最佳松弛因子。但好处是不用确定任何参数，他是对称正定线性方程组的方法也是求解大型稀疏线性方程组最热门的方法

Bicg迭代法本质是用双共轭梯度求解线性方程组的方法，对求解的方程没有正定性要求

Bicgstab迭代法本质是用稳定双共轭梯度求解线性方程组的方法，对求解的方程没有正定性要求

Gauss消去法第一、二个函数ltri、utri是一定要掌握的，后面的几乎每个函数都要用到 ltri 简单来说，当Ly=b b,L(非奇异下三角矩阵）已知求y

function y = ltri(L,b) n=size(b,1); y=zeros(n,1); for j = 1:n-1 y(j) = b(j)/L(j,j); b(j+1:n) = b(j+1:n) - y(j) * L(j+1:n,j); end y(n) = b(n)/L(n,n);

utri 简单来说，当Ux=y y,U(非奇异上三角矩阵）已知求x

function x = utri(U,y) n=size(y,1); x=zeros(n,1); for j = n:-1:2 x(j) = y(j)/U(j,j); y(1:j-1) = y(1:j-1) - x(j) * U(1:j-1,j); end x(1) = y(1)/U(1,1);

gauss算法，计算时粘贴过去就好

function [L,U] = gauss(A) n=size(A,1); for k = 1:n-1 A(k+1:n,k) = A(k+1:n,k) / A(k,k); A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)* A(k,k+1:n); end L = tril(A,-1) + eye(n); U = triu(A);

使用例子已经知道一个线性方程组，这里我就不写出数学形式了，A是系数矩阵，直接把上面写好的函数复制过来在运算就可以。主要是把这两个矩阵写进去。 A = [1 -1 1 -3;0 -1 -1 1;2 -2 -4 6;1 -2 -4 1]; b = [1;0;-1;-1];

A = [1 -1 1 -3;0 -1 -1 1;2 -2 -4 6;1 -2 -4 1]; b = [1;0;-1;-1]; [L,U] = gauss(A); y = ltri(L,b); x = utri(U,y)

Cholesky分解法这里是要用到ltri和utri，直接把上面的粘贴过来就好这里的算法比较简单，用L = chol(A,‘lower’);这里A是系数矩阵

使用例子用Cholesky分解法求接对称正定方程Ax=b，b随机选取，系数矩阵A是100阶的矩阵 n = 100; A = 10 * eye(n) + diag(ones(n-1,1),-1) + diag(ones(n-1,1),1);

n = 100; A = 10 * eye(n) + diag(ones(n-1,1),-1) + diag(ones(n-1,1),1); b = rand(n,1); L = chol(A,'lower'); y = ltri(L,b); x = utri(L',y)

Jacobi迭代法这个就是Jacobi迭代法的算法，用时复制。这里一共要输入四个条件，系数矩阵A，矩阵b，x0=zeros(列数,1),tol=1.0e-6

function x = jacobi(A,b,x0,tol) D = diag(diag(A)); L = -tril(A,-1); U = -triu(A,1); B = D \ (L + U); g = D \ b; x = B * x0 + g; n = 1; while norm(x - x0) > tol x0 = x; x= B * x0 + g; n = n + 1; end x,n

使用例子用Jacobi迭代法求解下列方程 A = [10 -1 0;-1 10 -2;0 -2 10]; b = [9;7;6];

A = [10 -1 0;-1 10 -2;0 -2 10]; b = [9;7;6]; x = zeros(3,1); tol = 1.0e-6; x = jacobi(A,b,x,tol);

Gauss-Seidel迭代法算法函数，输入四个条件，和上面相同

function x = seidel(A,b,x0,tol) D=diag(diag(A)); L = -tril(A,-1); U = -triu(A,1); B = (D - L) \ U; g = (D - L) \ b; x = B * x0 +g; n=1; while norm(x - x0) > tol x0 = x; x= B * x0 + g; n = n + 1; end x,n

用这个方法解上面Jacobi迭代法的例子

A = [10 -1 0;-1 10 -2;0 -2 10]; b = [9;7;6]; x = zeros(3,1); tol = 1.0e-6; x = seidel(A,b,x,tol);

超松弛迭代法算法函数，这里需要输入五个条件：系数矩阵A，矩阵b，x0=zeros(列数,1),w(松弛因子一到二之间）=1.025，tol=1.0e-6

function x = sor(A,b,x0,w,tol) D=diag(diag(A)); L = -tril(A,-1); U = -triu(A,1); B = (D - w * L) \ ((1 - w) * D + w * U); g = (D - w * L) \ b * w; x = B * x0 +g; n=1; while norm(x - x0) >= tol x0 = x; x= B * x0 + g; n = n + 1; end x,n

使用例子求这个方程组 A = [10 -1 0;-1 10 -2;0 -2 10]; b = [9;7;6];

A = [10 -1 0;-1 10 -2;0 -2 10]; b = [9;7;6]; x = zeros(3,1); tol = 1.0e-6; w = 1.025; x = sor(A,b,x,w,tol);

共轭梯度法

A = [16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12]; b = [32;26;38;30]; tol = 1.0e-6; maxit = 1000; x = pcg(A,b,tol,maxit)

Bicg迭代法这个比较简单，不需要额外写函数，直接用x = bicg(A,b,tol,maxit)maxit取1000

使用例子

n = 100; A = 10 * eye(n) + diag(ones(n-1,1),-1) + diag(ones(n-1,1),1); b = rand(n,1); tol = 1.0e-6; x = bicg(A,b,tol,1000)

Bicgstab迭代法和上面一样的输入条件 x = bicgstab(A,b,tol,1000)

n = 100; A = 10 * eye(n) + diag(ones(n-1,1),-1) + diag(ones(n-1,1),1); b = rand(n,1); tol = 1.0e-6; x = bicgstab(A,b,tol,1000)

【本文地址】

公司简介

联系我们