第三章 解线性方程组的直接方法的matlab程序

在科技、工程、医学和经济等各个领域中，经常遇到求解包含n 个未知数、由m 个方程构成的线性方程组

??????

?=+++=+++=+++m

n mn m m n n n n b x a x a x a b x a x a x a b x a x a x a ΛΛΛΛΛ22112

222212********* （3.1） 的问题.线性方程组（3.1）的解法一般有直接法和迭代法.迭代法是求解大型稀疏矩

阵方程组，尤其是由微分方程离散后得到的大型方程组的重要方法.在这章中我们要学习线性方程组的直接法，特别是适合用数学软件在计算机上求解的方法.

3.1 方程组的逆矩阵解法及其MATLAB 程序

3.1.1 逆矩阵、行列式及其MATLAB 命令

当矩阵A 为方阵时，A 的行列式通常表示为)det(A 或A .当0)det(≠A 时，A 可逆，且可以用表 3-1中的MATLAB 命令求A 的逆矩阵和行列式.

3.1.2 方程组的逆矩阵解法及其MATLAB 程序

在MATLAB 中引进了矩阵除法的概念，它包括左除和右除.设两矩阵B A ,为n 阶方阵，且A ，B 皆可逆，则矩阵方程C AX =，D YB =,F AZB =的解都可以用表 3-2中的命令求出.

3.1.3 线性方程组有解的判定条件及其MATLAB 程序

解非齐次线性方程组，将其增广矩阵化成行阶梯形矩阵，便可判断其是否有解．若

function [RA,RB,n]=jiepb(A,b)

B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0,

disp('请注意：因为RA~=RB ，所以此方程组无解.') return end

if RA==RB if RA==n

disp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') else

disp('请注意：因为RA=RB

例3.1.4 判断下列线性方程组解的情况.如果有唯一解，则用表 3-2方法求解.

(1) ???????=-+-=+-+=-++=+-+;

0742,

0634,

0723,

05324321432143214321x x x x x x x x x x x x x x x x (2) ??????

?=++-=+-+=-+-=+-+;0327,01613114,02332,075434321432143214321x x x x x x x x x x x x x x x x (3) ?????=+=+-=-+;8311,1023,22421

321321x x x x x x x x (4) ?????=--+=+-+=+-+.

12,2224,12w z y x w z y x w z y x 解 在MATLAB 工作窗口输入程序

>> A=[2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7]; b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b)

运行后输出结果为

请注意：因为RA=RB=n ，所以此方程组有唯一解. RA = 4,RB =4,n =4 在MATLAB 工作窗口输入

>>X=A\b,

运行后输出结果为 X =(0 0 0 0)’.

(2) 在MATLAB 工作窗口输入程序

>> A=[3 4 -5 7;2 -3 3 -2;4 11 -13 16;7 -2 1 3];b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b)

运行后输出结果

请注意：因为RA=RB

(3) 在MATLAB 工作窗口输入程序

>> A=[4 2 -1;3 -1 2;11 3 0]; b=[2;10;8]; [RA,RB,n]=jiepb(A,B)

运行后输出结果

请注意：因为RA~=RB ，所以此方程组无解. RA =2,RB =3,n =3

(4)在MATLAB 工作窗口输入程序

>> A=[2 1 -1 1;4 2 -2 1;2 1 -1 -1];

b=[1; 2; 1]; [RA,RB,n]=jiepb(A,b)

运行后输出结果

请注意：因为RA=RB

RA =2,RB =2,n =3

3.2 三角形方程组的解法及其MATLAB程序

用初等行变换解线性方程组是解线性方程组的一般方法.从例3.1.6可见，当方程的个数和未知量的个数很大时，用初等行变换解线性方程组的工作量相当的大，如果在计算过程中某一步发生错误，那么就会得出错误的结果.人们非常想用计算机的程序解线性方程组，而高斯（Gauss）消元法以及各种变形的高斯消元法使得用初等行变换解线性方程组实现机械化计算成为可能.

3.2.2 解三角形方程组的MATLAB程序

现提供的MATLAB程序是利用回代法解上三角形线性方程组，解下三角形线性方程

B=[A b]; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA;

if zhica>0,

disp('请注意：因为RA~=RB，所以此方程组无解.')

return

end

if RA==RB

if RA==n

disp('请注意：因为RA=RB=n，所以此方程组有唯一解.')

X=zeros(n,1); X(n)=b(n)/A(n,n);

for k=n-1:-1:1

X(k)=(b(k)-sum(A(k,k+1:n)*X(k+1:n)))/A(k,k);

end

else

disp('请注意：因为RA=RB

end

end

例3.2.2用解上三角形线性方程组的MATLAB程序解例3.2.1中的方程组.

解在MATLAB工作窗口输入程序

>>A=[5 -1 2 3;0 -2 7 -4;0 0 6 5;0 0 0 3];

b=[20; -7; 4;6];

[RA,RB,n,X]=shangsan(A,b)

运行后输出结果

请注意：因为RA=RB=n，所以此方程组有唯一解.

RA = RB =

4, 4,

n =

4,

X =[2.4 -4.0 -1.0 2.0]’

3.3 高斯（Gauss）消元法和列主元消元法及其MATLAB程序

3.3.1 高斯消元法及其MATLAB 程序

B=[A b]; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0,

disp('请注意：因为RA~=RB ，所以此方程组无解.') return end

if RA==RB if RA==n

disp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1

for k=p+1:n

m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);

end end

b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1

X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q); end

else

disp('请注意：因为RA=RB

例3.3.2 用高斯消元法和MATLAB 程序求解下面的非齐次线性方程组，并且用逆矩阵解方程组的方法验证.

??????

?-=+---=+--=+--=-+-.

142,16422,0,13432143214324321x x x x x x x x x x x x x x x 解 在MATLAB 工作窗口输入程序

>> A=[1 -1 1 -3; 0 -1 -1 1;2 -2 -4 6;1 -2 -4 1]; b=[1;0; -1;-1]; [RA,RB,n,X] =gaus (A,b)

运行后输出结果

请注意：因为RA=RB=n ，所以此方程组有唯一解. RA =

4

RB =

4

n =

4

3.3.2 列主元消元法及其MATLAB 程序

X = 0 -0.5000 0.5000 0

现提供的MATLAB程序是利用回代法解上三角形线性方程组.关于解下三角形线性

B=[A b]; n=length(b); RA=rank(A);

RB=rank(B);zhica=RB-RA;

if zhica>0,

disp('请注意：因为RA~=RB，所以此方程组无解.')

return

end

if RA==RB

if RA==n

disp('请注意：因为RA=RB=n，所以此方程组有唯一解.')

X=zeros(n,1); C=zeros(1,n+1);

for p= 1:n-1

[Y,j]=max(abs(B(p:n,p))); C=B(p,:);

B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;

for k=p+1:n

m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);

end

end

b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n);

for q=n-1:-1:1

X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);

end

else

disp('请注意：因为RA=RB

end

例3.3.3 用列主元消元法解线性方程组的MATLAB程序解例3.3.1的方程组.

解在MATLAB工作窗口输入程序

>> A=[0 -1 -1 1;1 -1 1 -3;2 -2 -4 6;1 -2 -4 1];

b=[0;1;-1;-1]; [RA,RB,n,X]=liezhu(A,b)

运行后输出结果

请注意：因为RA=RB=n，所以此方程组有唯一解.

RA = 4，RB = 4，n = 4，X =[0 -0.5 0.5 0]’

3.4 LU分解法及其MATLAB程序

3.4.1判断矩阵LU分解的充要条件及其MATLAB程序

[n n] =size(A); RA=rank(A);

if RA~=n

disp('请注意：因为A 的n 阶行列式hl 等于零，所以A 不能进行LU 分解.A 的秩

RA 如下：'), RA,hl=det(A); return

end

if RA==n

for p=1:n,h(p)=det(A(1:p, 1:p));, end

hl=h(1:n); for i=1:n

if h(1,i)==0

disp('请注意：因为A 的r 阶主子式等于零，所以A 不能进行LU 分解.A

的秩RA 和各阶顺序主子式值hl 依次如下：'),hl;RA,return

end end

if h(1,i)~=0

disp('请注意：因为A 的各阶主子式都不等于零，所以A 能进行LU 分解.A

的秩RA 和各阶顺序主子式值hl 依次如下：')

hl;RA end end

例3.4.1 判断下列矩阵能否进行LU 分解，并求矩阵的秩.

（1）????? ??6547121321；（2）????? ??654721321；（3）??

??? ??654321321. 解 （1）在MATLAB 工作窗口输入程序

>> A=[1 2 3;1 12 7;4 5 6];hl=pdLUfj(A)

运行后输出结果为

请注意：因为A 的各阶主子式都不等于零，所以A 能进行LU 分解.A 的秩RA 和

各阶顺序主子式值hl 依次如下：

RA = 3， hl = 1 10 -48

（2）在MATLAB 工作窗口输入程序

>> A=[1 2 3;1 2 7;4 5 6];hl=pdLUfj(A)

运行后输出结果为

请注意：因为A 的r 阶主子式等于零，所以A 不能进行LU 分解.A 的秩RA 和各阶

顺序主子式值hl 依次如下：

RA = 3， hl =1 0 12

（3）在MATLAB 工作窗口输入程序

>> A=[1 2 3;1 2 3;4 5 6];hl=pdLUfj(A)

运行后输出结果为

请注意：因为A 的n 阶行列式hl 等于零，所以A 不能进行LU 分解.A 的秩RA 如下

RA = 2， hl = 0

3.4.2 直接LU 分解法及其MATLAB 程序

[n n] =size(A); RA=rank(A); if RA~=n

disp('请注意：因为A 的n 阶行列式hl 等于零，所以A 不能进行LU 分解.A 的

秩RA 如下'), RA,hl=det(A);

return

end

if RA==n for p=1:n

h(p)=det(A(1:p, 1:p)); end

hl=h(1:n); for i=1:n

if h(1,i)==0

disp('请注意：因为A 的r 阶主子式等于零，所以A 不能进行LU 分解.A 的

秩RA 和各阶顺序主子式值hl 依次如下：'), hl;RA

return end end

if h(1,i)~=0

disp('请注意：因为A 的各阶主子式都不等于零，所以A 能进行LU 分解.A 的

秩RA 和各阶顺序主子式值hl 依次如下：')

for j=1:n

U(1,j)=A(1,j); end

for k=2:n for i=2:n for j=2:n

L(1,1)=1;L(i,i)=1; if i>j

L(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1); L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k))/U(k,k); else

U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j); end end end end

hl;RA,U,L end end

例3.4.3 用矩阵进行直接LU 分解的MA TLAB 程序分解矩阵=A ??????

?

?

?301

034

2

11010

0201.

解 在MATLAB 工作窗口输入程序

>> A=[1 0 2 0;0 1 0 1;1 2 4 3;0 1 0 3]; hl=zhjLU(A)

运行后输出结果

请注意：因为A 的各阶主子式都不等于零，所以A 能进行LU 分解.A 的秩RA

和各阶顺序主子式值hl 依次如下：

RA = 4

U = 1 0 2 0

0 1 0 1

0 0 2 1

0 0 0 2 由此可见，例3.4.3和例3.4.2中A 的LU 分解式相同.

3.4.3 含交换矩阵的LU 分解法及其MATLAB 程序

含交换矩阵的LU 分解法的MATLAB 程序是使用列主元消元法进行的，其命令和功能如下：

L = 1 0 0 0 0 1 0 0 1 2 1 0 0 1 0 1 hl = 1 1 2 4

3.4.4 判断正定对称矩阵的方法及其MATLAB 程序

[n n] =size(A); for p=1:n

h(p)=det(A(1:p, 1:p)); end

hl=h(1:n);zA=A'; for i=1:n

if h(1,i)<=0

disp('请注意：因为A 的各阶顺序主子式hl 不全大于零，所以A 不是正定

的.A 的转置矩阵zA 和各阶顺序主子式值hl 依次如下：'), hl;zA,return

end end

if h(1,i)>0

disp('请注意：因为A 的各阶顺序主子式hl 都大于零，所以A 是正定的.A 的转置

矩阵zA 和各阶顺序主子式值hl 依次如下：')

hl;zA end

例3.4.5 判断下列矩阵是否是正定对称矩阵：

（1）??

?

??

?? ?

?--9875411321114321

43

2

1.0；(2) ??

?

??

?? ?

?------196316902303112

11

； (3) ?

??

???

??

????

?

?--

--212100212100002121002121；（4）????

?

??---40106111

2. 解 （1）在MATLAB 工作窗口输入程序

>> A=[0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9];hl=zddc (A)

运行后输出结果

请注意：因为A 的各阶顺序主子式hl 不全大于零，所以A 不是正定的.A 的转

置矩阵zA 和各阶顺序主子式值hl 依次如下：

zA = 1/10 -1 11 5 2 2 21 7 3 -3 13 8 4 4 41 9 hl = 1/10 11/5 -1601/10 3696/5

因此,A 即不是正定矩阵，也不是对称矩阵.

（2）在MATLAB 工作窗口输入程序

>> A=[1 -1 2 1;-1 3 0 -3;2 0 9 -6;1 -3 -6 19],hl=zddc(A)

运行后输出结果

A = 1 -1 2 1

-1 3 0 -3

2 0 9 -6

1 -3 -6 19

请注意：因为A的各阶顺序主子式hl都大于零，所以A是正定的.A的转置矩阵zA 和各阶顺序主子式值hl依次如下：

zA = 1 -1 2 1

-1 3 0 -3

2 0 9 -6

1 -3 -6 19

hl = 1 2 6 24

（3）在MATLAB工作窗口输入程序

>> A=[1/sqrt(2) -1/sqrt(2) 0 0; -1/sqrt(2) 1/sqrt(2) 0 0; 0 0 1/sqrt(2) -1/sqrt(2); 0 0 -1/sqrt(2) 1/sqrt(2)], hl=zddc(A)

运行后输出结果

A= 985/1393 -985/1393 0 0

-985/1393 985/1393 0 0

0 0 985/1393 -985/1393

0 0 -985/1393 985/1393

请注意：因为A的各阶顺序主子式hl不全大于零，所以A不是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：

zA = 985/1393 -985/1393 0 0

-985/1393 985/1393 0 0

0 0 985/1393 -985/1393

0 0 -985/1393 985/1393

hl = 985/1393 0 0 0

可见，A不是正定矩阵，是半正定矩阵；因为A= A T因此，A是对称矩阵.

（4）在MATLAB工作窗口输入程序

>> A=[-2 1 1;1 -6 0;1 0 -4];hl=zddc(A)

运行后输出结果

A = -2 1 1

1 -6 0

1 0 -4

请注意：因为A的各阶顺序主子式hl不全大于零，所以A不是正定的.A的转置矩阵zA和各阶顺序主子式值hl依次如下：

zA = -2 1 1 hl = -2 11 -38

1 -6 0

1 0 -4

A不是正定矩阵，是负定矩阵；因为A= A T因此，A是对称矩阵.

3.4.5 正定对称矩阵的楚列斯基（Cholesky）分解及其MATLAB程序

>>help chol.

3.5 求解线性方程组的LU方法及其MATLAB程序

3.5.1 解线性方程组的楚列斯基（Cholesky ）分解法及其MATLAB 程序

例3.5.1 先将矩阵A 进行楚列斯基分解,然后解矩阵方程b AX =，并用其他方法验证.

??????

?

??=???

???? ??------=7531,19631690230311211b A . 解 根据（3.26）式编写MATLAB 程序，然后在工作窗口输入

>>A=[1 -1 2 1;-1 3 0 -3; 2 0 9 -6;1 -3 -6 19];

b1=1:2:7; b=b1'; R=chol(A);C=A-R'*R,R1=inv(R);R2=R1'; x=R1*R2*b,Rx=A\b-x

运行后输出方程组的解和验证结果

x = Rx = 1.0e-014 * C = 1.0e-015 *

-8.0000 -0.7105 0 0 0 0 0.3333 -0.0833 0 -0.4441 0 0 3.6667 0.2220 0 0 0 0 2.0000 0.1332 0 0 0 0

3.5.2 解线性方程组的直接LU 分解法及其MATLAB 程序 例3.5.2 首先将矩阵A 直接进行LU 分解，然后解矩阵方程b AX =

?????

??

?

?=30

1

0342110100201A ，??????

?

??-=5121b . 解 (1) 首先将矩阵A 直接进行LU 分解.在MATLAB 工作窗口输入程序

>> A=[1 0 2 0;0 1 0 1;1 2 4 3;0 1 0 3];b=[1;2;-1;5]; hl=zhjLU(A),A-L*U

运行后输出LU 分解

请注意：因为A 的各阶主子式都不等于零，所以A 能进行LU 分解.A 的秩RA

和各阶顺序主子式值hl 依次如下：

RA = 4

U = 1 0 2 0

0 1 0 1

0 0 2 1

0 0 0 2

A 分解为一个单位下三角形矩阵L 和一个上三角形矩阵

U 的积 LU A =. (2) 然后根据（3.27）式编写MATLAB 程序，然后在工作窗口输入

>> U=[1 0 2 0;0 1 0 1;0 0 2 1;0 0 0 2]; L=[1 0 0 0;0 1 0 0;1 2 1 0;0 1 0 1];

b=[1;2;-1;5];U1=inv(U); L1=inv(L); X=U1*L1*b,x=A\b

运行后输出方程组的解

X = 8.50000000000000 0.50000000000000 -3.75000000000000 1.50000000000000

3.5.3 解线性方程组的选主元的LU 方法及其MATLAB 程序

例3.5.3 先将矩阵A 进行LU 分解，然后解矩阵方程b AX = 其中

L = 1 0 0 0 0 1 0 0 1 2 1 0 0 1 0 1 hl = 1 1 2 4

???

???? ??--=98754113211143214321.0A ，??????

?

??-=5121b . 解 方法1 根据（3.28）式编写MATLAB 程序，然后在工作窗口输入

>> A=[0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9];

b=[1;2;-1;5]; [L U P]=LU(A), U1=inv(U); L1=inv(L); X=U1* L1*P*b

运行后输出结果

L = 1.0000 0 0 0 -0.0909 1.0000 0 0 0.0091 0.4628 1.0000 0 0.4545 -0.6512 0.2436 1.0000 U =11.0000 21.0000 13.0000 41.0000 0 3.9091 -1.8182 7.7273

0 0 3.7233 0.0512

0 0 0 -4.6171

方法2 根据（3.29）式编写MATLAB 程序，然后在工作窗口输入

>> A=[0.1 2 3 4;-1 2 -3 4;11 21 13 41;5 7 8 9];

b=[1;2;-1;5]; [F U]=LU(A), U1=inv(U); F1=inv(F); X=U1*F1*b

运行后输出结果

F=0.0091 0.4628 1.0000 0 -0.0909 1.0000 0 0 1.0000 0 0 0 0.4545 -0.6512 0.2436 1.0000 X =[-1.2013 3.3677 0.0536 -1.4440]’

解线性方程组A

b X =的选主元的LU 方法，现提供MATLAB 程序如下：

[n n] =size(A);B=[A b]; RA=rank(A); RB=rank(B); for p=1:n

h(p)=det(A(1:p, 1:p)); end

hl=h(1:n); for i=1:n

if h(1,i)==0

disp('请注意：因为A 的r 阶主子式等于零，所以A 不能进行LU 分解.A 的

秩RA 和各阶顺序主子式值hl 依次如下：')

hl;RA return end end

if h(1,i)~=0

disp('请注意：因为A 的各阶主子式都不等于零，所以A 能进行LU 分解.A

的秩RA 和各阶顺序主子式值hl 依次如下：')

X=zeros(n,1); Y=zeros(n,1); C=zeros(1,n);r=1:1; for p=1:n-1

[max1,j]=max(abs(A(p:n,p))); C=A(p,:); A(p,:)= A(j+P1,:); C= A(j+P1,:); g=r(p); r(p)= r(j+P1); r(j+P1)=g; for k=p+1:n

P = 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 1 X =[-1.2013 3.3677 0.0536 -1.4440]’ U=11.0000 21.0000 13.0000 41.0000 0 3.9091 -1.8182 7.7273 0 0 3.7233 0.0512 0 0 0 -4.6171

H= A(k,p)/A(p,p); A(k,p) = H; A(k,p+1:n)=A(k,p+1:n)- H* A(p,p+1:n); end end

Y(1)=B(r(1)); for k=2:n

Y(k)= B(r(k))- A(k,1:k-1)* Y(1:k-1); end

X(n)= Y(n)/ A(n,n); for i=n-1:-1:1

X(i)= (Y(i)- A(i, i+1:n) * X (i+1:n))/ A(i,i); end end

[RA,RB,n,X,Y]’;

3.6 误差分析及其两种MATLAB 程序

对于实际问题导出的线性方程组 b AX = 的系数矩阵A 和右端向量b 往往带有误差（扰动），下面讨论当A 或b 的微小变化对解X 的影响.

3.6.1 用MATLAB 软件作误差分析

例3.6.2 解下列矩阵方程b AX =，并比较方程（1）和（2）有何区别，它们的解有何变化.其中

,13/112/111/110/19/18/17

/112/111/110/19/18/17/16/111

/110

/19/18/17/16/15/110/19/18/17/16/15/14

/19/18/17/16/15/14/13/18/17/16/15/14/13/12/17/16

/15/14/13/12/11)

1(??

?

?

??

?

?

?

??

??=A ;

2222221?

?????????

? ??=b ,13/112/111/110/19/18/17/112/111/110/19/18/17/16/111/110/19/18/17/16/15/110/19/18/17/16/15/14/19/18/17/16/15/14/13/18/17/16/15/14/13/12/17/16/15/14/13/12/1001.1)

2(?

?

?

???

?

?

?

?

?

??=A .

2222221??????????? ??=b 解 (1) 矩阵方程b AX =的系数矩阵A 为7阶希尔伯特（Hilbert ）矩阵，我们可以用下列命令计算n 阶希尔伯特矩阵

>>h=hilb(n) % 输出h 为n 阶Hilbert 矩阵 在MATLAB 工作窗口输入程序

>> A=hilb(7);b=[1;2;2;2;2;2;2];X=A\b 运行后输出b AX =的解为 X =（-35 504 -1260 -4200 20790 -27720 12012T

).

（2）在MATLAB 工作窗口输入程序

>> B =[0.001,zeros(1,6);zeros(6,1),zeros(6,6)]; A=(B+hilb(7)); b=[1;2;2;2;2;2;2];X=A\b

运行后输出方程的解为 X =（-33 465 -966 -5181 22409 -29015 12413T

).

在MATLAB 工作窗口输入程序

>> X =[-33, 465,-966,-5181,22409,-29015,12413]';

X1 =[-35,504,-1260,-4200,20790,-27720,12012]'; wu=X1'- X'

运行后输出方程（1）和（2）的解的误差为

=

δX 401- 1295 1619- 981 294- 39 -2(T ).

方程（1）和（2）的系数矩阵的差为????

??=δ???661

661001.0O O O A ，常数向量相同，则b Ax =的解的差为=δX 4011295

1619981294392(----T ).A 的微小变化，

引起X 的很大变化，即X 对A 的扰动是敏感的.

3.6.2 求P 条件数和讨论b AX =解的性态的MATLAB 程序

根据定理3.10可知，判断线性方程组（3.38）是否病态需要计算条件数，在第一章中我们已经给出了计算矩阵A 的2范数条件数、1范数条件数、∞范数条件数、弗罗贝尼乌斯（Frobenius ）范数条件数、条件数倒数等MATLAB 命令.用这些命令计算各种条件数

Acw = cond (A, inf);Ac1= cond (A,1);

Ac2= cond (A,2);Acf = cond (A,'fro');dA=det(A); if (Acw>50)&(dA<0.1)

disp('请注意：AX=b 是病态的，A 的∞条件数,1条件数,2条件数,Frobenius 条件数和A 的行列式的值依次如下：') Acp=[Acw Ac1 Ac2 Acf dA]'; else

disp(' AX=b 是良态的，A 的∞条件数,1条件数,2条件数,Frobenius 条件数和A 的行列式的值依次如下：') Acp=[Acw Ac1 Ac2 Acf dA]'; end

n 阶希尔伯特矩阵定义为=H ?????

?

? ??-++)12/(1)1/(1/1)1/(13

/12/1/12/11n n n n n ΛM

M M ΛΛ，当n 很大时呈严重病态，常作为研究病态现象的系数矩阵.

例3.6.3 根据定理3.10，讨论线性方程组b AX =解的性态，并且求出A 的4种条件数.其中

（1）A 为7阶希尔伯特矩阵；（2）??????

?

?

?-----=74

2

163147213

5132A . 解 （1）首先将求P 条件数和讨论b AX =解的性态的MATLAB 程序保存名为

zpjxpb.m 的M 文件，然后在MATLAB 工作窗口输入程序

>> Acp =zpjxpb(hilb(7)); Acp',det(hilb(7))

运行后输出结果

请注意：AX=b 是病态的，A 的∞条件数,1条件数,2条件数,Frobenius 条件数

和A 的行列式的值依次如下：

ans = 1.0e+008 *

9.8519 9.8519 4.7537 4.8175 0.0000 ans = 4.8358e-025

（2）在MATLAB工作窗口输入程序

>> A=[2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7];Acp=zpjxpb(A); Acp'

运行后输出结果

AX=b是良态的，A的∞条件数,1条件数,2条件数,Frobenius条件数和A的行列式的值依次如下：

ans =

14.1713 19.4954 8.2085 11.4203 327.0000

3.6.3 用P范数讨论b

AX=解和A的性态的MATLAB程序

AX=（其中A是n阶希尔伯可以用下面提供的MATLAB程序讨论线性方程组b

AX=

特矩阵）解的性态，并利用（3.39）式和（3.40）式讨论当A或b有变化时，b

Acp = cond (A,P);dA=det(A); X=A\b;

dertaA=A-jA;

PndA=norm(dertaA, P);dertab=b-jb;Pndb=norm(dertab, P);

if Pndb>0

jX=A\jb; Pnb= norm(b, P);PnjX = norm(jX,P); dertaX=X-jX;

PnjdX= norm(dertaX, P);jxX= PnjdX/PnjX; PnjX = norm(jX,P);

PnX = norm(X,P); jxX= PnjdX/PnjX; xX= PnjdX/PnX;

Pndb=norm(dertab,P);

xAb=Pndb/Pnb;Pnbj=norm(jb,P); xAbj=Pndb/Pnbj;

Xgxx= Acp*xAb;

end

if PndA>0

jX=jA\b; dertaX=X-jX;PnX = norm(X,P);

PnjdX= norm(dertaX, P);

PnjX = norm(jX,P); jxX= PnjdX/PnjX;xX= PnjdX/PnX;

PnjA=norm(jA,P); PnA=norm(A,P);

PndA=norm(dertaA,P);xAbj= PndA/PnjA;xAb= PndA/PnA;

Xgxx= Acp*xAb;

end

if (Acp >50)&(dA<0.1)

disp('请注意：AX=b是病态的，A的P条件数Acp,A的行列式值dA，解X，近似解jX，解的相对误差jxX，解的相对误差估计值Xgxx，b或A的相对误差xAb依次如下：')

Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'

else

disp('请注意：AX=b是良态的，A的P条件数Acp,A的行列式值dA，解X，近似解jX，解的相对误差jxX，解的相对误差估计值Xgxx，b或A的相对误差xAb依次如下：')

Acp,dA,X',jX',xX',jxX',Xgxx',xAb',xAbj'

end

AX=解的性态，并利用（3.32）式例3.6.4 根据定理3.10，讨论线性方程组b

讨论当A的每个元都取4位有效数字时，其解的相对误差.其中A为7阶希尔伯特矩阵，

(

)

2

22243

1

1

=b T

.

解 （1）取∞范数和∞条件数，线性方程组b AX =的b 不变时，取∞范数和∞条

件数，系数矩阵A 为7阶希尔伯特矩阵，A 中的每个元素取4位有效数字，即

?

????

?????? ??=0.0769 0.0833 0.0909 0.1000 0.1111 0.1250 0.1429

0.0833 0.0909 0.1000 0.1111 0.1250 0.1429 0.1667 0.0909 0.1000 0.1111 0.1250 0.1429 0.1667 0.2000 0.1000 0.1111 0.1250 0.1429 0.1667 0.2000 0.2500 0.1111 0.1250 0.1429 0.1667 0.2000 0.2500 0.3333 0.1250 0.1429 0.1667 0.2000 0.2500 0.3333 0.5000 0.1429 0.1667 0.2000 0.2500 0.3333 0.5000 1.0000

jA . 用P 范数讨论b AX =解和A 的性态的MATLAB 程序保存名为zpjwc.m 的文件,然后

在工作窗口输入MATLAB 程序

>> jA =[1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.1429 0.1250 0.1111 0.1000 0.0909 0.0833 0.0769]; A=hilb(7);b=[1;1/3;4;2;2;2;2];

jb=[1;1/3;4;2;2;2;2]; Acp=zpjwc(A,jA,b,jb,inf)

运行后输出结果

请注意：AX=b 是病态的，A 的P 条件数Acp,A 的行列式值dA ，解X ，近似

解jX ，解的相对误差jxX ，解的相对误差估计值Xgxx ，b 或A 的相对误差xAb 依次如下：

Acp = dA =

9.8519e+008 4.8358e-025 ans =

1.0e+007 *

0.0020 -0.0697 0.6243 -2.3054 4.0677 -3.4123 1.0943 ans =

1.0e+004 *

0.0349 -0.4807 2.1126 -5.1087 7.6557 -6.3239 2.1112 xX = jxX = Xgxx =

0.9981 530.3248 4.9291e+004 xAb = xAbj = 5.0032e-005 5.0031e-005

由此可见，因为∞条件数(Cond ≈∞)A 985 194 889.719 848 31>>，所以此方程组为病态的b AX =的解

T 932)94210320,12334-790,676380,402.305-300,2436 256,697- 565,19( =X ，

b X jA =)(的解为

T

)11221239,63-,55767087,51-126,21 079,48- ,493( =jX 解的相对误差120.998≈X X δ,291.27

49≤+∞

∞

X X X δδ , 530.32≈+∞

∞X X X δδ, ,1025.0035-∞

∞

?≈A

A δ即相对误差放大了约985 194 889.72倍.

(2) 如果取2范数和2条件数计算，在MATLAB 工作窗口输入程序

>> Acp =zpjwc(A,jA,b,jb,2)

运行后输出结果

请注意：AX=b 是病态的，A 的P 条件数Acp,A 的行列式值dA ，解X ，近似

解jX ，解的相对误差jxX ，解的相对误差估计值Xgxx ，b 或A 的相对误差xAb 依次如下：

Acp = dA =

4.7537e+008 4.8358e-025 ans =

1.0e+007 *

0.0020 -0.0697 0.6243 -2.3054 4.0677 -3.4123 1.0943 ans =

1.0e+004 *

0.0349 -0.4807 2.1126 -5.1087 7.6557 -6.3239 2.1112 xX = jxX = Xgxx=

0.9981 511.0640 2.9951e+004 xAb = xAbj = 6.3006e-005 6.3005e-005

因为2条件数Cond ≈2)(A 475 367 356.591>>，所以此方程组为病态的.解的相对误差 10.998 22≈X X δ,,105300.6522-?≈A A δ,06.51122

≈+X X X δδ.85.9502922≤+X X X δδ 即相对误差放大了约475 367 356.59倍.

MatLab求解线性方程组

MatLab解线性方程组一文通 当齐次线性方程AX=0,rank(A)=r

MATLAB代码 解线性方程组的迭代法

解线性方程组的迭代法 1.rs里查森迭代法求线性方程组Ax=b的解 function[x,n]=rs(A,b,x0,eps,M) if(nargin==3) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值elseif(nargin==4) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-A)*x0+b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 2.crs里查森参数迭代法求线性方程组Ax=b的解 function[x,n]=crs(A,b,x0,w,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-w*A)*x0+w*b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x;

if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 3.grs里查森迭代法求线性方程组Ax=b的解 function[x,n]=grs(A,b,x0,W,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1;%前后两次迭代结果误差 %迭代过程 while(tol>eps) x=(I-W*A)*x0+W*b;%迭代公式 n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 4.jacobi雅可比迭代法求线性方程组Ax=b的解 function[x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3 eps=1.0e-6; M=200; elseif nargin<3 error return elseif nargin==5 M=varargin{1}; end D=diag(diag(A));%求A的对角矩阵 L=-tril(A,-1);%求A的下三角阵

matlab解方程组

matlab解方程组 lnx表示成log(x) 而lgx表示成log10(x) 1-exp(((log(y))/x^0.5)/(x-1)) 1、解方程 最近有多人问如何用matlab解方程组的问题，其实在matlab中解方程组还是很方便的，例如，对于代数方程组Ax=b(A为系数矩阵，非奇异)的求解，MATLAB 中有两种方法： (1)x=inv(A)*b —采用求逆运算解方程组； (2)x=A\B —采用左除运算解方程组 PS：使用左除的运算效率要比求逆矩阵的效率高很多~ 例: x1+2x2=8 2x1+3x2=13 >>A=[1,2;2,3];b=[8;13]; >>x=inv(A)*b x = 2.00 3.00 >>x=A\B x = 2.00 3.00； 即二元一次方程组的解x1和x2分别是2和3。 对于同学问到的用matlab解多次的方程组，有符号解法，方法是：先解出符号解,然后用vpa(F,n)求出n位有效数字的数值解.具体步骤如下： 第一步:定义变量syms x y z ...; 第二步:求解[x,y,z,...]=solve('eqn1','eqn2',...,'eqnN','var1','var2',...'varN'); 第三步:求出n位有效数字的数值解x=vpa(x,n);y=vpa(y,n);z=vpa(z,n);...。 如：解二（多）元二（高）次方程组： x^2+3*y+1=0 y^2+4*x+1=0 解法如下： >>syms x y; >>[x,y]=solve('x^2+3*y+1=0','y^2+4*x+1=0'); >>x=vpa(x,4); >>y=vpa(y,4); 结果是：

MATLAB解线性方程组的直接方法

在这章中我们要学习线性方程组的直接法，特别是适合用数学软件在计算机上求解的方法. 3.1 方程组的逆矩阵解法及其MATLAB 程序 3.1.3 线性方程组有解的判定条件及其MATLAB 程序 判定线性方程组A n m ?b X =是否有解的MATLAB 程序 function [RA,RB,n]=jiepb(A,b) B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB ，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') else disp('请注意：因为RA=RB> A=[2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7]; b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b) 运行后输出结果为 请注意：因为RA=RB=n ，所以此方程组有唯一解. RA = 4,RB =4,n =4 在MATLAB 工作窗口输入 >>X=A\b, 运行后输出结果为 X =(0 0 0 0)’. (2) 在MATLAB 工作窗口输入程序 >> A=[3 4 -5 7;2 -3 3 -2;4 11 -13 16;7 -2 1 3];b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b)

利用MATLAB求线性方程组

《MATLAB语言》课成论文 利用MATLAB求线性方程组 姓名：郭亚兰 学号：12010245331 专业：通信工程 班级：2010级通信工程一班 指导老师：汤全武 学院：物电学院 完成日期：2011年12月17日

利用MATLAB求解线性方程组 （郭亚兰 12010245331 2010 级通信一班） 【摘要】在高等数学及线性代数中涉及许多的数值问题，未知数的求解，微积分，不定积分，线性方程组的求解等对其手工求解都是比较复杂，而MATLAB语言正是处理线性方程组的求解的很好工具。线性代数是数学的一个分支，它的研究对象是向量，向量空间（或称线性空间），线性变换和有限维的线性方程组。因而，线性代数被广泛地应用于抽象代数和泛函分析中；由于科学研究中的非线性模型通常可以被近似为线性模型，使得线性代数被广泛地应用于自然科学和社会科学中。线性代数是数学的一个分支，它的研究对象是向量，向量空间（或称线性空间），线性变换和有限维的线性方程组。因而，线性代数被广泛地应用于抽象代数和泛函分析中；由于科学研究中的非线性模型通常可以被近似为线性模型，使得线性代数被广泛地应用于自然科学和社会科学中。线性代数是讨论矩阵理论、与矩阵结合的有限维向量空间及其线性变换理论的一门学科。 【关键字】线性代数MATLAB语言秩矩阵解 一、基本概念 1、N级行列式A：A等于所有取自不同性不同列的n个元素的积的代数和。 2、矩阵B：矩阵的概念是很直观的，可以说是一张表。 3、线性无关：一向量组（a1,a2,…,an）不线性相关，既没有不全为零的数 k1,k2,………kn使得：k1*a1+k2*a2+………+kn*an=0 4、秩：向量组的极在线性无关组所含向量的个数成为这个向量组的秩。 5、矩阵B的秩：行秩，指矩阵的行向量组的秩；列秩类似。记：R(B)

线性方程组求解matlab实现

3.1 方程组的逆矩阵解法及其MATLAB 程序 3.1.3 线性方程组有解的判定条件及其MATLAB 程序 判定线性方程组A n m ?b X =是否有解的MATLAB 程序 function [RA,RB,n]=jiepb(A,b) B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB ，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') else disp('请注意：因为RA=RB> A=[2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7]; b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b) 运行后输出结果为 请注意：因为RA=RB=n ，所以此方程组有唯一解. RA = 4,RB =4,n =4 在MATLAB 工作窗口输入 >>X=A\b, 运行后输出结果为 X =(0 0 0 0)’. (2) 在MATLAB 工作窗口输入程序 >> A=[3 4 -5 7;2 -3 3 -2;4 11 -13 16;7 -2 1 3];b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b) 运行后输出结果 请注意：因为RA=RB> A=[4 2 -1;3 -1 2;11 3 0]; b=[2;10;8]; [RA,RB,n]=jiepb(A,B) 运行后输出结果 请注意：因为RA~=RB ，所以此方程组无解. RA =2,RB =3,n =3 (4)在MATLAB 工作窗口输入程序

实验一用matlab求解线性方程组

实验1.1 用matlab 求解线性方程组 第一节 线性方程组的求解 一、齐次方程组的求解 rref （A ） %将矩阵A 化为阶梯形的最简式 null （A ） %求满足AX ＝0的解空间的一组基，即齐次线性方程组的基 础解系 【例１】 求下列齐次线性方程组的一个基础解系，并写出通解： 我们可以通过两种方法来解： 解法1： >> A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; >> rref(A) 执行后可得结果： ans= 1 -1 0 0 0 0 -1 1 0 0 0 0 由最简行阶梯型矩阵，得化简后的方程 ??? ??=+--=+--=-+-0 22004321 43214321x x x x x x x x x x x x

取x2，x4为自由未知量，扩充方程组为 即 提取自由未知量系数形成的列向量为基础解系，记 所以齐次方程组的通解为 解法2： clear A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; B=null(A, 'r') % help null 看看加个‘r’是什么作用, 若去掉r ，是什么结果？ 执行后可得结果： B= 1 0 1 0 0 1 0 1 ?? ?=-=-0 04321x x x x ?????? ?====4 4432221x x x x x x x x ??? ??? ??????+????????????=????? ???????1100001142 4321x x x x x x , 00111????? ? ??????=ε, 11002????? ???????=ε2 211εεk k x +=

Matlab线性方程组求解(Gauss消去法)

Matlab线性方程组求解 1. Gauss消元法： function x=DelGauss(a,b) % Gauss消去法 [n,m]=size(a); nb=length(b); det=1; %存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if a(k,k)==0 return end 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); 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 Example: >> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]'; >> x=DelGauss(A,b) x = 0.9739 -0.0047 1.0010 2. 列主元Gauss消去法： function x=detGauss(a,b) % Gauss列主元消去法 [n,m]=size(a); nb=length(b); det=1; %存储行列式值 x=zeros(n,1); for k=1:n-1 amax=0; %选主元 for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k));r=i; end end if amax<1e-10 return; end if r>k %交换两行 for j=k:n

Matlab求解线性方程组非线性方程组

求解线性方程组 solve，linsolve 例： A=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1]; %矩阵的行之间用分号隔开，元素之间用逗号或空格 B=[3;1;1;0] X=zeros(4,1);%建立一个4元列向量 X=linsolve(A,B) diff（fun，var，n）：对表达式fun中的变量var求n阶导数。 例如：F=sym（'u(x,y)*v(x,y)'）; %sym（）用来定义一个符号表达式 diff(F); %matlab区分大小写 pretty(ans) %pretty（）：用习惯书写方式显示变量；ans是答案表达式 非线性方程求解 fsolve(fun,x0,options) 为待解方程或方程组的文件名；fun其中 x0位求解方程的初始向量或矩阵； option为设置命令参数 建立文件fun.m： function y=fun(x) y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ... x(2) - 0.5*cos(x(1))+0.3*sin(x(2))]; >>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve')) 注： ...为续行符 m文件必须以function为文件头，调用符为@；文件名必须与定义的函数名相同；fsolve（）主要求解复杂非线性方程和方程组，求解过程是一个逼近过程。Matlab求解线性方程组 AX=B或XA=B 在MATLAB中，求解线性方程组时，主要采用前面章节介绍的除法运算符“/”和“\”。如： X=A\B表示求矩阵方程AX＝B的解； 的解。XA=B表示矩阵方程B/A＝X． 对方程组X＝A\B，要求A和B用相同的行数，X和B有相同的列数，它的行数等于矩阵A的列数，方程X＝B/A同理。 如果矩阵A不是方阵，其维数是m×n，则有： m＝n 恰定方程，求解精确解； m>n 超定方程，寻求最小二乘解； m

线性方程组求解Matlab程序(精.选)

线性方程组求解 1.直接法 Gauss消元法： function x=DelGauss(a,b) % Gauss消去法 [n,m]=size(a); nb=length(b); det=1;%存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if a(k,k)==0 return end 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); 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 Example: >> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]'; >> x=DelGauss(A,b) x = 0.9739 -0.0047 1.0010 列主元Gauss消去法： function x=detGauss(a,b) % Gauss列主元消去法

[n,m]=size(a); nb=length(b); det=1;%存储行列式值 x=zeros(n,1); for k=1:n-1 amax=0;% 选主元 for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k));r=i; end end if amax<1e-10 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

matlab常用解方程及方程组函数

1. roots 求解多项式的根 r=roots(c) 注意： c 为一维向量，者返回指定多项式的所有根( 包括复根),poly 和roots 是互为反运算，还有就是roots 只能求解多项式的解 还有下面几个函数poly2sym、sym2poly 、eig >>syms x >>y=x A5+3*x A3+3; >>c=sym2poly(y);%求解多项式系数 >>r=roots(c); >>poly(r) 2. residue 求留数 [r, p, k] = residue(b,a) >>b = [ 5 3 -2 7] >>a = [-4 0 8 3] >>[r, p, k] = residue(b,a) 3. solve 符号解方程(组)——使用最多的 g = solve(eq1,eq2,...,eqn,var1,var2,...,varn) 注意：eqn 和varn 可以是符号表达式，也可以是字符串表达式，但是使用符号表达式时不能有“=号”，假如说varn 没有给出，使用findsym 函数找出默认的求解变量。返回的g 是个结构体，以varn 为字段。由于符号求解的局限性，好多情况下可能得到空矩阵，此时只能用数值解法 解方程A=solve('a*xA2 + b*x + c') 解方程组B=solve('a*uA2 + vA2', 'u - v = 1', 'aA2 - 5*a + 6') 4. fzero 数值求零点 [x,fval,exitflag,output]=fzero(fun,x0,options,p1,p2...) fun 是目标函数，可以是句柄(@)、inline 函数或M 文件名 x0 是初值，可以是标量也可以是长度为2 的向量，前者给定一个位置，后者是给定一个范围options 是优化参数，通过optimset 设置，optimget 获取，一般使用默认的就可以了，具体参照帮助 p1,p2...为需要传递的其它参数 假如说(x/1446)A2+p/504.1+(t/330.9)*(log(1-x/1446)+(1-1 /5.3)*x/1446)=0 的根，其中p,t 是已知

用matlab解线性方程组

用matlab解线性方程组 2008-04-12 17:00 一。高斯消去法 1.顺序高斯消去法 直接编写命令文件 a=[] d=[]' [n,n]=size(a); c=n+1 a(:,c)=d; for k=1:n-1 a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去 end x=[0 0 0 0]' %回带 x(n)=a(n,c)/a(n,n); for g=n-1:-1:1 x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g) end 2.列主高斯消去法 *由于“[r,m]=max(abs(a(k:n,k)))”返回的行是“k:n,k”内的第几行，所以要通过修正来把m 改成真正的行的值。该程序只是演示程序，真正机器计算不需要算主元素所在列以下各行应为零的值。 直接编写命令文件 a=[] d=[] ' [n,n]=size(a); c=n+1 a(:,c)=d; %(增广) for k=1:n-1 [r,m]=max(abs(a(k:n,k))); %选主 m=m+k-1; %(修正操作行的值) if(a(m,k)~=0) if(m~=k) a([k m],:)=a([m k],:); %换行 end a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k))*a(k, k:c); %消去end end x=[0 0 0 0]' %回带 x(n)=a(n,c)/a(n,n); for g=n-1:-1:1 x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g) end

雅可比解线性方程组matlab

雅可比迭代 使用雅可比迭代法求解线性方程组的步骤 步骤1：输入系数矩阵A和方程组右端向量B； 步骤2：将矩阵A分解为下三角阵L对角阵D和上三角阵U 可分解为（D+L+U）X=B for o=1:n d(o,o)=a(o,o); u(o,o+1:n)=-a(o,o+1:n); end for p=2:n l(p,1:p-1)=-a(p,1:p-1); end； 步骤3：将上式化简为x=B0x+f,其中B0=-D-1(L+U),f=D-1B for i=1:n b0(i,i+1:n)=u(i,i+1:n)/a(i,i); f(i,:)=b(i,:)/a(i,i); end for p=2:n b0(p,1:p-1)=l(p,1:p-1)/a(p,p);；

步骤4：采用迭代公式在允许误差范围e=1e-7内求得解向量x x0=x; x=b0*x+f 雅可比迭代法matlab程序： function [x,k]=jacobi(a,b) n=length(a); e=1e-7; m=100; x0=zeros(n,1); x=x0; k=0; d=zeros(n); l=zeros(n); u=zeros(n); b0=zeros(n); f=zeros(n,1);

x0=x+2*e; for o=1:n d(o,o)=a(o,o); u(o,o+1:n)=-a(o,o+1:n); end for p=2:n l(p,1:p-1)=-a(p,1:p-1); end for i=1:n b0(i,i+1:n)=u(i,i+1:n)/a(i,i); f(i,:)=b(i,:)/a(i,i); end for p=2:n b0(p,1:p-1)=l(p,1:p-1)/a(p,p); end while max(abs(x0-x))>e&k

解线性方程组直接方法matlab用法

2.1 方程组地逆矩阵解法及其MATLAB 程序 2.1.3 线性方程组有解地判定条件及其MATLAB 程序判定线性方程组A n m ?b X =是否有解地MATLAB 程序 function [RA,RB,n]=jiepb(A,b) B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB ，所以此方程 组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n ，所以此方 程组有唯一解.') else disp('请注意：因为RA=RB

Matlab求解线性方程组、非线性方程组

Matlab求解线性方程组、非线性方程组 姓名：罗宝晶学号：15 专业：材料学院高分子系 第一部分数值计算 Matlab求解线性方程组AX=B或XA=B 在MATLAB中，求解线性方程组时，主要采用除法运算符“/”和“\”。如：X=A\B表示求矩阵方程AX＝B的解； X＝B/A表示矩阵方程XA=B的解。 对方程组X＝A\B，要求A和B用相同的行数，X和B有相同的列数，它的行数等于矩阵A的列数，方程X＝B/A同理。 如果矩阵A不是方阵，其维数是m×n，则有： m＝n 恰定方程，求解精确解； m>n 超定方程，寻求最小二乘解； mm。则方程组没有精确解，此时称方程组为超定方程组。线性超定方程组经常遇到的问题是数据的曲线拟合。对于超定方程，在MATLAB中，利用左除命令（x=A\b）来寻求它的最小二乘解；还可以用广义逆来求，即x=pinv(A),所得的解不一定满足Ax=b，x只是最小二乘意义上的解。左除的方法是建立在奇异值分解基础之上，由此获得的解最可靠；广义逆法是建立在对原超定方程直接进行householder变换的基础上，其算法可靠性稍逊与奇异值求解，但速度较快；

matlab常用解方程及方程组函数

matlab常用解方程及方程组函数 1.roots求解多项式的根 r=roots(c) 注意：c为一维向量，者返回指定多项式的所有根(包括复根),poly和roots是互为反运算，还有就是roots只能求解多项式的解 还有下面几个函数poly2sym、sym2poly、eig >>syms x >>y=x^5+3*x^3+3; >>c=sym2poly(y);%求解多项式系数 >>r=roots(c); >>poly(r) 2.residue求留数 [r, p, k] = residue(b,a) >>b = [ 5 3 -2 7] >>a = [-4 0 8 3] >>[r, p, k] = residue(b,a) 3.solve符号解方程(组)——使用最多的 g = solve(eq1,eq2,...,eqn,var1,var2,...,varn) 注意：eqn和varn可以是符号表达式，也可以是字符串表达式，但是使用符号表达式时不能有“=”号，假如说varn没有给出，使用findsym函数找出默认的求解变量。返回的g是一个结构体，以varn为字段。由于符号求解的局限性，好多情况下可能得到空矩阵，此时只能用数值解法 解方程A=solve('a*x^2 + b*x + c') 解方程组B=solve('a*u^2 + v^2', 'u - v = 1', 'a^2 - 5*a + 6') 4.fzero数值求零点 [x,fval,exitflag,output]=fzero(fun,x0,options,p1,p2...) fun是目标函数，可以是句柄(@)、inline函数或M文件名 x0是初值，可以是标量也可以是长度为2的向量，前者给定一个位置，后者是给定一个范围 options是优化参数，通过optimset设置，optimget获取，一般使用默认的就可以了，具体参照帮助 p1,p2...为需要传递的其它参数 假如说(x/1446)^2+p/504.1+(t/330.9)*(log(1-x/1446)+(1-1/5.3)*x/1446)=0的根，其中p,t是已知参数，但是每次都改变 那么目标函数如下三种书写格式，效果完全等效。注意参数列表中，未知数一定放第一位，其他参数放后面 (1)objfun=@(x,p,t)(x/1446).^2+p/504.1+(t/330.9).*(log(1-x/1446)+(1-1/5.3).*x/1446); (2)objfun=inline('(x/1446).^2+p/504.1+(t/330.9).*(log(1-x/1446)+(1-1/5.3).*x/1446)','x','p','t') 此时的调用格式如下 fzero(objfun,x0,options,p,t)%如果options使用的默认的话，那直接使用[]，p和t就是我们需要传递的参数 fzero(@(x)objfun(x,p,t),x0,options)%这种格式与上面的等效 区别就是前者，将参数p和t作为fzero的参数进行传递，而后者是将p和t作为objfun的参数进行传递，没有本质区别 (3)function f=objfun(x,p,t)%以M文件格式书写目标函数 f=(x/1446).^2+p/504.1+(t/330.9).*(log(1-x/1446)+(1-1/5.3).*x/1446); 此时有三种调用格式

MATLAB之GAUSS消元法解线性方程组

Matlab之Gauss消元法解线性方程组 1.Gauss消元法 function x=DelGauss(a,b) %Gauss消去法 [n,m]=size(a); nb=length(b); det=1;%存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if a(k,k)==0 return end 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); 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 Example: >>A=[1.0170-0.00920.0095;-0.00920.99030.0136;0.00950.0136 0.9898]; >>b=[101]'; >>x=DelGauss(A,b) x= 0.9739 -0.0047 1.0010 2.列主元Gauss消去法： function x=detGauss(a,b) %Gauss列主元消去法 [n,m]=size(a); nb=length(b); det=1;%存储行列式值 x=zeros(n,1); for k=1:n-1 amax=0;%选主元 for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k));r=i; end end if amax<1e-10 return;

Matlab求解线性方程组、非线性方程组.docx

求解线性方程组solve ，linsolve 例： A=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1]; %矩阵的行之间用分号隔开，元素之间用逗号或空格 B=[3;1;1;0] X=zeros(4,1);% 建立一个4 元列向量 X=linsolve(A,B) diff (fun , Var, n)：对表达式fun中的变量Var求n阶导数。 例如：F=sym('u(x,y)*v(x,y)' ) ; %sym ()用来定义一个符号表达式diff(F); %matlab 区分大小写 pretty(ans) %pretty ()：用习惯书写方式显示变量；ans 是答案表达式非线性方程求解 fsolVe(fun,x0,options) 其中fun 为待解方程或方程组的文件名； x0 位求解方程的初始向量或矩阵； option 为设置命令参数 建立文件fun.m ： function y=fun(x) y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ... x(2) - 0.5*cos(x(1))+0.3*sin(x(2))]; >>clear;x0=[0.1,0.1];fsolVe(@fun,x0,optimset('fsolVe')) 注： ...为续行符 m 文件必须以function 为文件头,调用符为@；文件名必须与定义的函数名相同；fsolVe ()主要求解复杂非线性方程和方程组,求解过程是一个逼近过程。 Matlab 求解线性方程组 AX=B 或XA=B

在MATLAB 中，求解线性方程组时，主要采用前面章节介绍的除法运算符“和/ ” “”。如： X=A?B表示求矩阵方程AX = B的解； X= B/A表示矩阵方程XA=B的解。 对方程组X = A?B ,要求A和B用相同的行数，X和B有相同的列数，它的行数等于矩阵A 的列数，方程X= B/A 同理。 如果矩阵A不是方阵，其维数是m× n,则有：m = n 恰定方程，求解精确解；m>n 超定方程，寻求最小二乘解； m

实验一用matlab求解线性方程组

实验 用matlab 求解线性方程组 第一节 线性方程组的求解 一、齐次方程组的求解 rref （A ） %将矩阵A 化为阶梯形的最简式 null （A ） %求满足AX ＝0的解空间的一组基，即齐次线性方程组的基 础解系 【例１】 求下列齐次线性方程组的一个基础解系，并写出通解： 我们可以通过两种方法来解： 解法1： >> A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; >> rref(A) 执行后可得结果： ans= 1 -1 0 0 0 0 -1 1 0 0 0 0 由最简行阶梯型矩阵，得化简后的方程 ??? ??=+--=+--=-+-0 22004321 43214321x x x x x x x x x x x x ?? ?=-=-0 04321x x x x

取x2，x4为自由未知量，扩充方程组为 即 提取自由未知量系数形成的列向量为基础解系，记 所以齐次方程组的通解为 解法2： clear A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; B=null(A, 'r') % help null 看看加个‘r’是什么作用,若去掉 r ，是什么结果 执行后可得结果： B= 1 0 1 0 0 1 0 1 易见，可直接得基础解系 ?????? ?====4 4432221x x x x x x x x ??? ??? ??????+????????????=????? ???????1100001142 4321x x x x x x , 00111????? ? ??????=ε, 11002????? ???????=ε2 211εεk k x +=, 0111? ?????????=ε, 1002? ?????????=ε

线性代数方程组数值解法及MATLAB实现综述

线性代数方程组数值解法及MATLAB 实现综述 廖淑芳 20122090 数计学院 12计算机科学与技术1班（职教本科） 一、分析课题 随着科学技术的发展，提出了大量复杂的数值计算问题，在建立电子计算机成为数值计算的主要工具以后，它以数字计算机求解数学问题的理论和方法为研究对象。其数值计算中线性代数方程的求解问题就广泛应用于各种工程技术方面。因此在各种数据处理中，线性代数方程组的求解是最常见的问题之一。关于线性代数方程组的数值解法一般分为两大类：直接法和迭代法。 直接法就是经过有限步算术运算，可求的线性方程组精确解的方法（若计算过程没有舍入误差），但实际犹如舍入误差的存在和影响，这种方法也只能求得近似解，这类方法是解低阶稠密矩阵方程组级某些大型稀疏矩阵方程组的有效方法。直接法包括高斯消元法，矩阵三角分解法、追赶法、平方根法。 迭代法就是利用某种极限过程去逐步逼近线性方程组精确解的方法。迭代法具有需要计算机的存储单元少，程序设计简单，原始系数矩阵在计算过程始终不变等优点，但存在收敛性级收敛速度问题。迭代法是解大型稀疏矩阵方程组（尤其是微分方程离散后得到的大型方程组）的重要方法。迭代法包括Jacobi 法SOR 法、SSOR 法等多种方法。 二、研究课题-线性代数方程组数值解法 一、 直接法 1、 Gauss 消元法 通过一系列的加减消元运算，也就是代数中的加减消去法，以使A 对角线以下的元素化为零，将方程组化为上三角矩阵；然后，再逐一回代求解出x 向量。 1.1消元过程 1. 高斯消元法(加减消元)：首先将A 化为上三角阵，再回代求解。 11121121222212n n n n nn n a a a b a a a b a a a b ?? ? ? ? ???L L M M O M M L (1)(1)(1)(1)(1)11121311(2)(2)(2)(2)222322(3)(3)(3)3333()()000000n n n n n nn n a a a a b a a a b a a b a b ?? ? ? ? ? ? ???L L L M M M O M M L 步骤如下：