心得Word文档下载推荐.docx

心得Word文档下载推荐.docx

《心得Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《心得Word文档下载推荐.docx（44页珍藏版）》请在冰豆网上搜索。

Digit_Matrix=[1/3sqrt

（2）3.4234；

exp（0.23）log（29）23^（-11.23）]

Syms_Matrix=sym（Digit_Matrix）

结果是：

Digit_Matrix=

0.33331.41423.4234

1.25863.36730.0000

Syms_Matrix=

[1/3，sqrt

（2），17117/5000]

[5668230535726899*2^（-52），7582476122586655*2^（-51），5174709270083729*2^（-103）]

注意：

矩阵是用分数形式还是浮点形式表示的，将矩阵转化成符号矩阵后，都将以最接近原值的有理数形式表示或者是函数形式表示。

（4）符号简化

函数simple或simplify%寻找符号矩阵或符号表达式的最简型

格式simple（s）%s是矩阵或表达式

[R,how]=simple（s）%R为返回的最简形，how为简化过程中使用的主要方法。

说明Simple（s）将表达式s的长度化到最短。

若还想让表达式更加精美，可使用函数Pretty。

格式Pretty（s）%使表达式s更加精美

例1-64计算行列式

的值。

在Matlab编辑器中建立M文件：

symsabcd

A=[1111;

abcd;

a^2b^2c^2d^2;

a^4b^4c^4d^4];

d1=det（A）

d2=simple（d1）%化简表达式d1

pretty（d2）%让表达式d2符合人们的书写习惯

则显示结果如下：

d1=

b*c^2*d^4-b*d^2*c^4-b^2*c*d^4+b^2*d*c^4+b^4*c*d^2-b^4*d*c^2-a*c^2*d^4+a*d^2*c^4+a*b^2*d^4-a*b^2*c^4-a*b^4*d^2+a*b^4*c^2+a^2*c*d^4-a^2*d*c^4-a^2*b*d^4+a^2*b*c^4+a^2*b^4*d-a^2*b^4*c-a^4*c*d^2+a^4*d*c^2+a^4*b*d^2-a^4*b*c^2-a^4*b^2*d+a^4*b^2*c

d2=

（-d+c）*（b-d）*（b-c）*（-d+a）*（a-c）*（a-b）*（a+c+d+b）

（-d+c）（b-d）（b-c）（-d+a）（a-c）（a-b）（a+c+d+b）

1.4线性方程的组的求解

我们将线性方程的求解分为两类：

一类是方程组求唯一解或求特解，另一类是方程组求无穷解即通解。

可以通过系数矩阵的秩来判断：

若系数矩阵的秩r=n（n为方程组中未知变量的个数），则有唯一解；

若系数矩阵的秩r<

n，则可能有无穷解；

线性方程组的无穷解=对应齐次方程组的通解+非齐次方程组的一个特解；

其特解的求法属于解的第一类问题，通解部分属第二类问题。

1.4.1求线性方程组的唯一解或特解（第一类问题）

这类问题的求法分为两类：

一类主要用于解低阶稠密矩阵——直接法；

另一类是解大型稀疏矩阵——迭代法。

1．利用矩阵除法求线性方程组的特解（或一个解）

方程：

AX=b

解法：

X=A\b

例1-76求方程组

的解。

解：

A=[56000

15600

01560

00156

00015];

B=[10001]'

;

R_A=rank（A）%求秩

X=A\B%求解

运行后结果如下

R_A=

5

X=

2.2662

-1.7218

1.0571

-0.5940

0.3188

这就是方程组的解。

或用函数rref求解：

C=[A,B]%由系数矩阵和常数列构成增广矩阵C

R=rref（C）%将C化成行最简行

R=

1.000000002.2662

01.0000000-1.7218

001.0000001.0571

0001.00000-0.5940

00001.00000.3188

则R的最后一列元素就是所求之解。

例1-77求方程组

的一个特解。

A=[11-3-1;

3-1-34;

15-9-8];

B=[140]'

X=A\B%由于系数矩阵不满秩，该解法可能存在误差。

X=[00-0.53330.6000]’（一个特解近似值）。

若用rref求解，则比较精确：

A=[11-3-1;

C=[A,B];

%构成增广矩阵

R=rref（C）

1.00000-1.50000.75001.2500

01.0000-1.5000-1.7500-0.2500

00000

由此得解向量X=[1.2500–0.250000]’（一个特解）。

2．利用矩阵的LU、QR和cholesky分解求方程组的解

（1）LU分解：

LU分解又称Gauss消去分解，可把任意方阵分解为下三角矩阵的基本变换形式（行交换）和上三角矩阵的乘积。

即A=LU，L为下三角阵，U为上三角阵。

则：

A*X=b变成L*U*X=b

所以X=U\（L\b）这样可以大大提高运算速度。

命令[L，U]=lu（A）

例1-78求方程组

A=[42-1;

3-12;

1130];

B=[2108]'

D=det（A）

[L,U]=lu（A）

X=U\（L\B）

显示结果如下：

D=

0

L=

0.3636-0.50001.0000

0.27271.00000

1.000000

U=

11.00003.00000

0-1.81822.0000

000.0000

Warning:

Matrixisclosetosingularorbadlyscaled.

Resultsmaybeinaccurate.RCOND=2.018587e-017.

InD:

\Matlab\pujun\lx0720.matline4

1.0e+016*

-0.4053

1.4862

1.3511

说明结果中的警告是由于系数行列式为零产生的。

可以通过A*X验证其正确性。

（2）Cholesky分解

若A为对称正定矩阵，则Cholesky分解可将矩阵A分解成上三角矩阵和其转置的乘积，即：

其中R为上三角阵。

方程A*X=b变成

所以

（3）QR分解

对于任何长方矩阵A，都可以进行QR分解，其中Q为正交矩阵，R为上三角矩阵的初等变换形式，即：

A=QR

方程A*X=b变形成QRX=b

所以X=R\（Q\b）

上例中[Q,R]=qr（A）

X=R\（Q\B）

说明这三种分解，在求解大型方程组时很有用。

其优点是运算速度快、可以节省磁盘空间、节省内存。

1.4.2求线性齐次方程组的通解

在Matlab中，函数null用来求解零空间，即满足A·

X=0的解空间，实际上是求出解空间的一组基（基础解系）。

格式z=null%z的列向量为方程组的正交规范基，满足

。

%z的列向量是方程AX=0的有理基

例1-79求解方程组的通解：

A=[1221;

21-2-2;

1-1-4-3];

formatrat%指定有理式格式输出

B=null（A,'

r'

）%求解空间的有理基

运行后显示结果如下：

B=

25/3

-2-4/3

10

01

或通过行最简行得到基：

B=rref（A）

1.00000-2.0000-1.6667

01.00002.00001.3333

0000

即可写出其基础解系（与上面结果一致）。

写出通解：

symsk1k2

X=k1*B（:

1）+k2*B（:

2）%写出方程组的通解

pretty（X）%让通解表达式更加精美

运行后结果如下：

[2*k1+5/3*k2]

[-2*k1-4/3*k2]

[k1]

[k2]

%下面是其简化形式

[2k1+5/3k2]

[]

[-2k1-4/3k2]

[k1]

[k2]

1.4.3求非齐次线性方程组的通解

非齐次线性方程组需要先判断方程组是否有解，若有解，再去求通解。

因此，步骤为：

第一步：

判断AX=b是否有解，若有解则进行第二步

第二步：

求AX=b的一个特解

第三步：

求AX=0的通解

第四步：

AX=b的通解=AX=0的通解+AX=b的一个特解。

例1-80求解方程组

在Matlab中建立M文件如下：

A=[1-23-1;

3-15-3;

212-2];

b=[123]'

B=[Ab];

n=4;

R_A=rank（A）

R_B=rank（B）

formatrat

ifR_A==R_B&

R_A==n%判断有唯一解

X=A\b

elseifR_A==R_B&

R_A<

n%判断有无穷解

X=A\b%求特解

C=null（A,'

）%求AX=0的基础解系

elseX='

equitionnosolve'

%判断无解

end

运行后结果显示：

2

R_B=

3

equitionnosolve

说明该方程组无解

例1-81求解方程组的通解：

解法一：

在Matlab编辑器中建立M文件如下：

b=[140]'

R_A==n

n

Equationhasnosolves'

运行后结果显示为：

2

Rankdeficient,rank=2tol=8.8373e-015.

\Matlab\pujun\lx0723.matline11

0

-8/15

3/5

C=

3/2-3/4

3/27/4

10

01

所以原方程组的通解为X=k1

+k2

+

解法二：

用rref求解

C=rref（B）%求增广矩阵的行最简形，可得最简同解方程组。

10-3/23/45/4

01-3/2-7/4-1/4

00000

对应齐次方程组的基础解系为：

，

非齐次方程组的特解为：

所以，原方程组的通解为：

X=k1

1.4.4线性方程组的LQ解法

函数symmlq

格式x=symmlq（A,b）%求线性方程组AX=b的解X。

A必须为n阶对称方阵，b为n元列向量。

A可以是由afun定义并返回A*X的函数。

如果收敛，将显示结果信息；

如果收敛失败，将给出警告信息并显示相对残差norm（b-A*x）/norm（b）和计算终止的迭代次数。

symmlq（A,b,tol）%指定误差tol，默认值是1e-6。

symmlq（A,b,tol,maxit）%maxit指定最大迭代次数

symmlq（A,b,tol,maxit,M）%M为用于对称正定矩阵的预处理因子

symmlq（A,b,tol,maxit,M1,M2）%M=M1×

M2

symmlq（A,b,tol,maxit,M1,M2,x0）%x0为初始估计值，默认值为0。

[x,flag]=symmlq（A,b,…）%flag的取值为：

0表示在指定迭代次数之内按要求精度收敛；

1表示在指定迭代次数内不收敛；

2表示M为坏条件的预处理因子；

3表示两次连续迭代完全相同；

4表示标量参数太小或太大；

5表示预处理因子不是对称正定的。

[x,flag,relres]=symmlq（A,b,…）%relres表示相对误差norm（b-A*x）/norm（b）

[x,flag,relres,iter]=symmlq（A,b,…）%iter表示计算x的迭代次数

[x,flag,relres,iter,resvec]=symmlq（A,b,…）%resvec表示每次迭代的残差：

norm（b-A*x0）

[x,flag,relres,iter,resvec,resveccg]=symmlq（A,b,…）%resveccg表示每次迭代共轭梯度残差的范数

1.4.5双共轭梯度法解方程组

函数bicg

格式x=bicg（A,b）%求线性方程组AX=b的解X。

A必须为n阶方阵，b为n元列向量。

bicg（A,b,tol）%指定误差tol，默认值是1e-6。

bicg（A,b,tol,maxit）%maxit指定最大迭代次数

bicg（A,b,tol,maxit,M）%M为用于对称正定矩阵的预处理因子

bicg（A,b,tol,maxit,M1,M2）%M=M1×

bicg（A,b,tol,maxit,M1,M2,x0）%x0为初始估计值，默认值为0。

[x,flag]=bicg（A,b,…）%flag的取值为：

4表示标量参数太小或太大。

[x,flag,relres]=bicg（A,b,…）%relres表示相对误差norm（b-A*x）/norm（b）

[x,flag,relres,iter]=bicg（A,b,…）%iter表示计算x的迭代次数

[x,flag,relres,iter,resvec]=bicg（A,b,…）%resvec表示每次迭代的残差：

例1-83调用MATLAB6.0数据文件west0479。

loadwest0479

A=west0479;

%将数据取为系数矩阵A。

b=sum（A,2）;

%将A的各行求和，构成一列向量。

X=A\b;

%用“\”求AX=b的解。

norm（b-A*X）/norm（b）%计算解的相对误差。

ans=

1.2454e-017

[x,flag,relres,iter,resvec]=bicg（A,b）%用bicg函数求解。

x=（全为0，由于太长，不显示出来）

flag=

1%表示在默认迭代次数（20次）内不收敛。

relres=%相对残差relres=norm（b-A*x）/norm（b）=norm（b）/norm（b）=1。

1

iter=%表明解法不当，使得初始估计值0向量比后来所有迭代值都好。

resvec=（略）%每次迭代的残差。

semilogy（0:

20,resvec/norm（b）,'

-o'

）%作每次迭代的相对残差图形，结果如下图。

xlabel（'

iterationnumber'

）%x轴为迭代次数。

ylabel（'

relativeresidual'

）%y轴为相对残差。

图1-1双共轭梯度法相对误差图

1.4.6稳定双共轭梯度方法解方程组

函数bicgstab

格式x=bicgstab（A,b）

bicgstab（A,b,tol）

bicgstab（A,b,tol,maxit）

bicgstab（A,b,tol,maxit,M）

bicgstab（A,b,tol,maxit,M1,M2）

bicgstab（A,b,tol,maxit,M1,M2,x0）

[x,flag]=bicgstab（A,b,…）

[x,flag,relres]=bicgstab（A,b,…）

[x,flag,relres,iter]=bicgstab（A,b,…）

[x,flag,relres,iter,resvec]=bicgstab（A,b,…）

稳定双共轭梯度法解方程组，调用方式和返回的结果形式和命令bicg一样。

例1-84

loadwest0479;

A=west0479;

b=sum（A,2）;

[x,flag]=bicgstab（A,b）

显示结果是x的值全为0，flag=1。

表示在默认误差和默认迭代次数（20次）下不收敛。

若改为：

[L1,U1]=luinc（A,1e-5）;

[x1,flag1]=bicgstab（A,b,1e-6,20,L1,U1）

即指定误差，并用A的不完全LU分解因子L和U作为预处理因子M=L*U，其结果是x1的值全为0，flag=2表示预处理因子为坏条件的预处理因子。

若改为

[L2,U2]=luinc（A,1e-6）;

%稀疏矩阵的不完全LU分解。

[x2,flag2,relres2,iter2,resvec2]=bicgstab（A,b,1e-15,10,L2,U2）

%指定最大迭代次数为10次，预处理因子M=L*U。

semilogy（0:

0.5:

iter2,resvec2/norm（b）,'

）%每次迭代的相对残差图形，见图1-2。

xlabel（'

ylabel（'

结果为

x2=（其值全为1，略）

flag2=

0%表示收敛

relres2=

2.8534e-016%收敛时的相对误差

iter2=

6%计算终止时的迭代次数

resvec2=%每次迭代的残差

1.4.7复共轭梯度平方法解方程组

函数cgs

格式x=cgs（A,b）

cgs（A,b,tol）

cgs（A,b,tol,maxit）

cgs（A,b,tol,maxit,M）

cgs（A,b,tol,maxit,M1,M2）

cgs（A,b,tol,maxit,M1,M2,x0）

[x,flag]=cgs（A,b,…）

[x,flag,relres]=cgs（A,b,…）

[x,flag,relres,iter]=cgs（A,b,…）

[x,flag,relres,iter,resvec]=cgs（A,b,…）

调用方式和返回的结果形式与命令bicg一样。

1.4.8共轭梯度的LSQR方法

函数lsqr

格式x=lsqr（A,b）

lsqr（A,b,tol）

lsqr（A,b,tol,maxit）

lsqr（A,b,tol,maxit,M）

lsqr（A,b,tol,maxit,M1,M2）

lsqr（A,b,tol,maxit,M1,M2,x0）

[x,flag]=lsqr（A,b,…）

[x,flag,relres]=lsqr（A,b,…）

[x,flag,relres,iter]=lsqr（A,b,…）

[x,flag,relres,iter,resvec]=lsqr（A,b,…）

例1-85

n=100;

on=ones（n,1）;

A=spdiags（[-2*on4*on-on],-1:

1,n,n）;

%产生一个对角矩阵

b=sum（A,2）;

tol=1e-8;

%指定精度

maxit=15;

%指定最大迭代次数

M1=spdiags（[on/（-2）on],-1:

0,n,n）;

M2=spdiags（[4*on-on],0:

%M1*M2=M，即产生预处理因子

[x,flag,relres,iter,resvec]=lsqr（A,b,tol,maxit,M1,M2,[]）

结果显示

x的值全为1

relres=

3.5241e-009%表示相对残差

iter=

12%计算终止时的迭代次数

1.4.9广义最小残差法

函数qmres

格式x=gmres（A,b）

gmres（A,b,restart）

gmres（A,b,restart,tol）

gm