MATLAB求解线性方程组

不会炼丹的炼丹术士

已于 2023-04-16 19:35:58 修改

阅读量5.1k

点赞数 8

文章标签： matlab 算法线性代数

于 2023-04-10 21:22:20 首次发布

本文链接：https://blog.csdn.net/m0_62490077/article/details/130054012

版权

1.使用左除或者是右除的方法

使用左除的话解决的是方程Ax=b的解

>> A=[3, 6, 5;4, 8, 3;2, 5, 1];
>> b=[7 5 9]';
>> x=A\b

x =

  -13.8182
    7.0909
    1.1818

使用右除的话，解决的是方程xB=D的解

>> B=[3, 6, 5;4, 8, 3;2, 5, 1];
>> D=[7 5 9];
>> x=D/B

x =

   -0.2727    6.4545   -9.0000

注意点1：

这里要注意计算时候各个矩阵的维度，比如在第一个代码块中， A是3*3的矩阵，x是3*1的矩阵，

b是3*1的矩阵，所以不要忘记丢掉转置符号，否则MATLAB会报错显示维度不对。

注意点2：

• 若方程组无解,除法给出最小二乘意义上的近似解，即使向量 AX －B的长度达到最小；

• 若方程组有无穷多解，除法将给出一个具有最多零元素的特解；

• 若为唯一解，除法将给出解。

2.使用求逆矩阵的方法

这里利用的是对于方程Ax=b，x=A^(-1)*b

对应的MATLAB函数是inv函数

>> A=[3, 6, 5;4, 8, 3;2, 5, 1];
>> b=[7 5 9]';
>> x=inv(A)*b

x =

  -13.8182
    7.0909
    1.1818

注意点1：

这种方法需要注意的是，只有方阵才有逆矩阵，否则会报错

>> A=[3, 6, 5;4, 8, 3;2, 5, 1;7,8,1];
>> b=[1,2,5,8]';
>> x=inv(A)*b
错误使用 inv
矩阵必须为方阵。

注意点2：

还有需要注意A不能是奇异的

>> A=[1,0,0;0,1,0;0,0,0];
>> b=[1,2,3]'
>> x=inv(A)*b
警告: 矩阵在工作精度内为奇异的。 
 

x =

     1
     2
   Inf

3. 用linsolve函数求解

>> A=[3 1 -1;1 2 4;-1 4 5];
b=[3.6;2.1;-1.4];
x=linsolve(A,b)

x =

    1.4818
   -0.4606
    0.3848

注意：

当 A 是方阵时， linsolve 使用 LU 分解和部分主元消元法；对于所有其他情况， linsolve 使用 QR 分解和列主元消元法。

4.RLS方法（递归最小二乘法）

迭代公式：

$x\left ( k \right )=x\left ( k-1 \right )+K\left ( k \right )\left ( y\left ( k \right )-b\left ( k \right ) ^\mathsf{T}x\left ( k-1 \right )\right )$

$K\left ( k \right )=\frac{P\left ( k-1\right )b\left ( k \right )}{1+b\left ( k \right )^\mathsf{T}P\left ( k-1\right )b\left ( k \right )}$

$P\left ( k \right )=\left [ I-K\left ( k \right )b\left ( k \right ) ^\mathsf{T}\right ]P\left ( k-1 \right )$

MATLAB实现代码如下

第一块是实现迭代部分的函数

function res=RLS(A,b)
%brief：求解Ax=b
%input  para：A,方程组Ax=b中A
%input  para： b，方程组Ax=b中b
%output para： res为辨识的参数，即求解得到的x值
format long;           %默认long型
[len,num] = size(A);   %len为A的行数(观测数据个数） ，num为A的列数
x = rand(num,1);       %设定初值
I = eye(num, num);     %单位矩阵，与观测数据维数相同
P = (10^6) * I;        %设定P的初始值，一般给定一个较大的值

%循环迭代部分
for k = 1:len
Ak = A(k,:);              %新的数据部分，A(k,:)表示的是矩阵A的第k行的所有元素，
                          %所以在这里指的是矩阵A的第k行
Q1 = P*(Ak');             %K(k)的分子部分
Q2 = 1 + Ak * P * (Ak');  %K(k)的分母部分
K = Q1/Q2;
x = x + K * (b(k) - Ak*x); %x的更新方程
P = (I - K*Ak)*P;          %P的更新方程

thetae(:,k) = x;           %A(:,k)表示矩阵A的第k列，这里thetae矩阵是用来保存每次迭代的x的值，每个x都是列向量
end                        %循环结束

res=(thetae(:,end));       %表示输出thetae矩阵的最后一列
end                        %函数结束

第二块是主函数部分

close all;
clear;
clc;
%3x1+6x2+5x3=7
%4x1+8x2+3x3=5
%2x1+5x2+1x3=9
%矩阵Ax=b, x=[x1;x2;x3]
%得A， b为
A=[3, 6, 5;4, 8, 3;2, 5, 1];
b=[7 5 9]';
ans = RLS(A,b)   %调用函数，计算结果

得出的解如下：

ans =

 -13.818053079745313
   7.090852901317834
   1.181806202286966

5.利用Kaczmarz算法求解

Kaczmarz特点是无需计算 $A^{T}A$ 的逆矩阵，该算法在矩阵A的行数（观测数据）非常多的情况下非常实用。

令 $a_{j}^{T}$ 表示矩阵A的第j行， $b_j$ 表示向量b的第j个元素，μ为正实数，满足0<μ<2.

Kaczmarz算法的步骤如下：

（1）令i=0，选定初始值 $x^{(0)}$ 。

（2）对于j=1，…，m，令

$x^{(im+j)}=x^{(im+j-1)}+\mu \left (b{_j}-a{_j}^{T}x^{(im+j-1)}\right )a{_j}/a{_j}^{T}a{_j}$

（3）令i=i+1，回到第二步。

下面是代码部分

第一块是函数的迭代函数的代码。

function [x,iter,error] = kaczmarz(A,b, x0,maxit,exactx,tol, miu)
% Ax = b 
% x0是初始值
% maxit最大迭代次数, tol为算法精度, exactx为真实的x值
% x - 所求的值
% iter - 迭代次数
% error - x和exactxd每次迭代时候的差值
% miu-每一步的步长
m = size(A,1);  % m表示A的行数
% n = size(A,2);n表示A的列数
x = x0;
error = [];
if isempty(tol) % 如果未给定误差精度，就按照给定的最大迭代次数进行
iter = maxit;
% 下面是循环部分
for i=1:maxit
pickedi = mod(i,m)+1;  % 返回m除以i后的余数再加1
row = A(pickedi, :);   % row表示的是A的第pickedi行，也就是aj'
x = x + miu*( b(pickedi) - (row * x) ) / (row * row') * row';
e = norm(x-exactx);
error = [error,e];
end
else
iter = 1; %如果给定误差精度，按照算法精度更新 
e = 1;
    while e >= tol
    pickedi = mod(iter,m)+1;
    row = A(pickedi, :);
    x = x + miu*( b(pickedi) - (row * x) ) / (row * row') * row';
    e = norm(x-exactx);
    error = [error,e];
    iter = iter+1;
    end
end
end

有时间更新随机化Kaczmarz算法。