直接解法
Ax=l，求x。利用逆矩阵、矩阵除法、x=inv(A'A)*A'l
lind.m文件如下，
%求解线性方程组
% x1 + 2x2 + 3x3 + 9x4 = 5
%2x1 + 2x2 + 5x3 + 4x4 = 2
%3x1 + 5x2 + x3 + 5x4 = 3
%7x1 + 4x2 + 2x3 -10x4 = 8
A=[1 2 3 9;2 2 5 4;3 5 1 5;7 4 2 -10];
b=[5 2 3 8]';
x=A\b
c=inv(A)*b
d=inv(A'*A)*A'*b
输出
>> lind.m
x =
4.9227
-3.2373
-1.1680
1.1173
c =
4.9227
-3.2373
-1.1680
1.1173
d =
4.9227
-3.2373
-1.1680
1.1173
Warning: Direct access of structure fields returned by a function call (e.g.,
call to lind) is not allowed. See MATLAB 7.10 Release Notes, "Subscripting Into Function Return Values" for details.
??? Attempt to reference field of non-structure array.
求线性方程组通解
lind.m内容如下
%求解线性方程组通解
% X1 +5 *X2 -1 *X3 -1 *X4 = -1
% X1 -2 *X2 +1 *X3 +3 *X4 = 3
%3*X1 +8 *X2 -1 *X3 +1 *X4 = 1
% X1 -9 *X2 +3 *X3 +7 *X4 = 7
%系数矩阵和增广矩阵的秩都为2，说明有多组解，且基础解系4-2=2
%求基础解系用Z=null(A,'r')，可获得标准化的基础解系。特解采用广义逆函数pinv(A)，
%可求得A的广义逆矩阵，采用x0=pinv(A)*b求得特解，通解为x=k1*x1+k2*x2+...+kn*xn+x0
A=[1 5 -1 -1 ;
1 -2 1 3;
3 8 -1 1;
1 -9 3 7];
b=[-1 3 1 7]';
r1=rank(A)%系数矩阵的秩
r2=rank(A,b)%增广矩阵的秩
Z=null(A,'r')%求基础解系
x0=pinv(A)*b%求得特解
%验证该解
k1=rand(1);
k2=rand(1);
x=k1*Z(:,1)+k2*Z(:,2)+x0
err=norm(A*x-b)
>> lind.m
r1 =
2
r2 =
2
Z =
-0.4286 -1.8571
0.2857 0.5714
1.0000 0
0 1.0000
x0 =
0.3785
-0.0876
0.1873
0.7530
x =
-0.1805
0.0994
0.2848
1.0315
err =
1.8971e-015
矛盾方程组最小二乘解
%矛盾方程组的最小二乘解
%原始数据表
% x 1 1 2 2 2 3 3 4 4 4 5 6
% y 14.8 15.9 20.2 20.0 18.55 22.2 20.9 21 18.3 20.7 16.1 14.75
% 分别一次
A=[ 1 1 2 2 2 3 3 4 4 4 5 6]';
t=ones(size(A));
A1=[t,A];
y= [14.8 15.9 20.2 20.0 18.55 22.2 20.9 21 18.3 20.7 16.1 14.75]';
%x=inv(A'*A)*A'*y
x1=pinv(A1)*y%线性回归方程系数
n=length(A);
rmse1=sqrt(sum((y-A1*x1).^2)/n)%均方误差
xsq=A.^2;
A2=[t,A,xsq]
x2=pinv(A2)*y%线性回归方程系数
rmse2=sqrt(sum((y-A2*x2).^2)/n)%均方误差
xmin=min(A);
xmax=max(A);
nx1=xmin:0.2:xmax;
ny1=x1(1)+x1(2)*nx1;
nx2=nx1;
ny2=x2(1)+x2(2)+x2(3)*nx2.^2;
plot(A',y,'*',nx1,ny1,'+b',nx2,ny2,'r')
输出
>> lind.m
x1 =
18.8820
-0.0861
rmse1 =
2.5105
A2 =
1 1 1
1 1 1
1 2 4
1 2 4
1 2 4
1 3 9
1 3 9
1 4 16
1 4 16
1 4 16
1 5 25
1 6 36
x2 =
10.3924
6.3994
-0.9793
rmse2 =
1.0960
