心得.docx

心得.docx

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

心得

1.1.2符号矩阵的生成

在MATLAB中输入符号向量或者矩阵的方法和输入数值类型的向量或者矩阵在形式上很相像，只不过要用到符号矩阵定义函数sym，或者是用到符号定义函数syms，先定义一些必要的符号变量，再像定义普通矩阵一样输入符号矩阵。

1．用命令sym定义矩阵：

这时的函数sym实际是在定义一个符号表达式，这时的符号矩阵中的元素可以是任何的符号或者是表达式，而且长度没有限制，只是将方括号置于用于创建符号表达式的单引号中。

如下例：

例1-3

>>sym_matrix=sym（'[abc；Jack，HelpMe!

，NOWAY!

]，'）

sym_matrix=

[abc]

[JackHelpMe!

NOWAY!

]

>>sym_digits=sym（'[123；abc；sin（x）cos（y）tan（z）]'）

sym_digits=

[123]

[abc]

[sin（x）cos（y）tan（z）]

2．用命令syms定义矩阵

先定义矩阵中的每一个元素为一个符号变量，而后像普通矩阵一样输入符号矩阵。

例1-4

>>symsabc；

>>M1=sym（'Classical'）；

>>M2=sym（'Jazz'）；

>>M3=sym（'Blues'）

>>syms_matrix=[abc；M1，M2，M3；int2str（[235]）]

syms_matrix=

[abc]

[ClassicalJazzBlues]

[235]

把数值矩阵转化成相应的符号矩阵。

数值型和符号型在MATLAB中是不相同的，它们之间不能直接进行转化。

MATLAB提供了一个将数值型转化成符号型的命令，即sym。

例1-5

>>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

线性方程组的无穷解=对应齐次方程组的通解+非齐次方程组的一个特解；其特解的求法属于解的第一类问题，通解部分属第二类问题。

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;3-1-34;15-9-8];

B=[140]';

>>C=[A,B];%构成增广矩阵

>>R=rref（C）

R=

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

X=

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）

B=

1.00000-2.0000-1.6667

01.00002.00001.3333

0000

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

写出通解：

symsk1k2

X=k1*B（:

1）+k2*B（:

2）%写出方程组的通解

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

运行后结果如下：

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

X=A\b%求特解

C=null（A,'r'）%求AX=0的基础解系

elseX='equitionnosolve'%判断无解

end

运行后结果显示：

R_A=

2

R_B=

3

X=

equitionnosolve

说明该方程组无解

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

解法一：

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

A=[11-3-1;3-1-34;15-9-8];

b=[140]';

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

X=A\b

C=null（A,'r'）

elseX='Equationhasnosolves'

end

运行后结果显示为：

R_A=

2

R_B=

2

Warning:

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

>InD:

\Matlab\pujun\lx0723.matline11

X=

0

0

-8/15

3/5

C=

3/2-3/4

3/27/4

10

01

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

+k2

+

解法二：

用rref求解

A=[11-3-1;3-1-34;15-9-8];

b=[140]';

B=[Ab];

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

运行后结果显示为：

C=

10-3/23/45/4

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

00000

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

，

非齐次方程组的特解为：

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

X=k1

+k2

+

。

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元列向量。

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

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

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×M2

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

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

0表示在指定迭代次数之内按要求精度收敛；1表示在指定迭代次数内不收敛；2表示M为坏条件的预处理因子；3表示两次连续迭代完全相同；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表示每次迭代的残差：

norm（b-A*x0）

例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向量比后来所有迭代值都好。

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）,'-o'）%每次迭代的相对残差图形，见图1-2。

>>xlabel（'iterationnumber'）

>>ylabel（'relativeresidual'）

结果为

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,…）

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

例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:

1,n,n）;%M1*M2=M，即产生预处理因子

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

结果显示

x的值全为1

flag=

0%表示收敛

relres=

3.5241e-009%表示相对残差

iter=

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

1.4.9广义最小残差法

函数qmres

格式x=gmres（A,b）

gmres（A,b,restart）

gmres（A,b,restart,tol）

gm