|
|
Matlab学习主题
MATLAB是美国MathWorks公司自20世纪80年代中期推出的数学软件，优秀的数值计算能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出。到目前为止，其最高版本6.1版已经推出。随着版本的不断升级，它在数值计算及符号计算功能上得到了进一步完善。MATLAB已经发展成为多学科、多种工作平台的功能强大的大型软件。在欧美等高校，MATLAB已经成为线性代数、自动控制理论、概率论及数理统计、数字信号处理、时间序列分析、动态系统仿真等高级课程的基本教学工具，是攻读学位的大学生、硕士生、博士生必须掌握的基本技能。
学一学吧,不管你是哪方面的研究可以都用到它.
——————————————————————————————————————————————
线性规划问题
线性规划问题是目标函数和约束条件均为线性函数的问题，MATLAB6.0解决的线性规划问题的标准形式为：
min
sub.to：
其中f、x、b、beq、lb、ub为向量，A、Aeq为矩阵。
其它形式的线性规划问题都可经过适当变换化为此标准形式。
在MATLAB6.0版中，线性规划问题（Linear Programming）已用函数linprog取代了MATLAB5.x版中的lp函数。当然，由于版本的向下兼容性，一般说来，低版本中的函数在6.0版中仍可使用。
函数 linprog
格式 x = linprog(f,A,b) %求min f ' *x sub.to 线性规划的最优解。
x = linprog(f,A,b,Aeq,beq) %等式约束 ，若没有不等式约束 ，则A=[ ]，b=[ ]。
x = linprog(f,A,b,Aeq,beq,lb,ub) %指定x的范围 ，若没有等式约束 ，则Aeq=[ ]，beq=[ ]
x = linprog(f,A,b,Aeq,beq,lb,ub,x0) %设置初值x0
x = linprog(f,A,b,Aeq,beq,lb,ub,x0,options) % options为指定的优化参数
[x,fval] = linprog(…) % 返回目标函数最优值，即fval= f ' *x。
[x,lambda,exitflag] = linprog(…) % lambda为解x的Lagrange乘子。
[x, lambda,fval,exitflag] = linprog(…) % exitflag为终止迭代的错误条件。
[x,fval, lambda,exitflag,output] = linprog(…) % output为关于优化的一些信息
说明 若exitflag>0表示函数收敛于解x，exitflag=0表示超过函数估值或迭代的最大数字，exitflag<0表示函数不收敛于解x；若lambda=lower 表示下界lb，lambda=upper表示上界ub，lambda=ineqlin表示不等式约束，lambda=eqlin表示等式约束，lambda中的非0元素表示对应的约束是有效约束；output=iterations表示迭代次数，output=algorithm表示使用的运算规则，output=cgiterations表示PCG迭代次数。
例5-1 求下面的优化问题
min
sub.to
解：
>>f = [-5; -4; -6];
>>A = [1 -1 1;3 2 4;3 2 0];
>>b = [20; 42; 30];
>>lb = zeros(3,1);
>>[x,fval,exitflag,output,lambda] = linprog(f,A,b,[],[],lb)
结果为：
x = %最优解
0.0000
15.0000
3.0000
fval = %最优值
-78.0000
exitflag = %收敛
1
output =
iterations: 6 %迭代次数
cgiterations: 0
algorithm: 'lipsol' %所使用规则
lambda =
ineqlin: [3x1 double]
eqlin: [0x1 double]
upper: [3x1 double]
lower: [3x1 double]
>> lambda.ineqlin
ans =
0.0000
1.5000
0.5000
>> lambda.lower
ans =
1.0000
0.0000
0.0000
表明：不等约束条件2和3以及第1个下界是有效的
————————————————————————————————————————————
非线性规划问题
有约束的一元函数的最小值
单变量函数求最小值的标准形式为 sub.to
在MATLAB5.x中使用fmin函数求其最小值。
函数 fminbnd
格式 x = fminbnd(fun,x1,x2) %返回自变量x在区间 上函数fun取最小值时x值，fun为目标函数的表达式字符串或MATLAB自定义函数的函数柄。
x = fminbnd(fun,x1,x2,options) % options为指定优化参数选项
[x,fval] = fminbnd(…) % fval为目标函数的最小值
[x,fval,exitflag] = fminbnd(…) %xitflag为终止迭代的条件
[x,fval,exitflag,output] = fminbnd(…) % output为优化信息
说明 若参数exitflag>0，表示函数收敛于x，若exitflag=0，表示超过函数估计值或迭代的最大数字，exitflag<0表示函数不收敛于x；若参数output=iterations表示迭代次数，output=funccount表示函数赋值次数，output=algorithm表示所使用的算法。
--------------------------------------------------------------------------------
线性方程的组的求解
1.4 线性方程的组的求解
我们将线性方程的求解分为两类：一类是方程组求唯一解或求特解，另一类是方程组求无穷解即通解。可以通过系数矩阵的秩来判断：
若系数矩阵的秩r=n（n为方程组中未知变量的个数），则有唯一解；
若系数矩阵的秩r<n，则可能有无穷解；
线性方程组的无穷解 = 对应齐次方程组的通解+非齐次方程组的一个特解；其特解的求法属于解的第一类问题，通解部分属第二类问题。
1.4.1 求线性方程组的唯一解或特解（第一类问题）
这类问题的求法分为两类：一类主要用于解低阶稠密矩阵 —— 直接法；另一类是解大型稀疏矩阵 —— 迭代法。
1．利用矩阵除法求线性方程组的特解（或一个解）
方程：AX=b
解法：X=A\\b
例1-76 求方程组 的解。
解：
>>A=[5 6 0 0 0
1 5 6 0 0
0 1 5 6 0
0 0 1 5 6
0 0 0 1 5];
B=[1 0 0 0 1]';
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.0000 0 0 0 0 2.2662
0 1.0000 0 0 0 -1.7218
0 0 1.0000 0 0 1.0571
0 0 0 1.0000 0 -0.5940
0 0 0 0 1.0000 0.3188
则R的最后一列元素就是所求之解。
例1-77 求方程组 的一个特解。
解：
>>A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
>>B=[1 4 0]';
>>X=A\\B %由于系数矩阵不满秩，该解法可能存在误差。
X =[ 0 0 -0.5333 0.6000]’（一个特解近似值）。
若用rref求解，则比较精确：
>> A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
B=[1 4 0]';
>> C=[A,B]; %构成增广矩阵
>> R=rref(C)
R =
1.0000 0 -1.5000 0.7500 1.2500
0 1.0000 -1.5000 -1.7500 -0.2500
0 0 0 0 0
由此得解向量X=[1.2500 – 0.2500 0 0]’（一个特解）。
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=[4 2 -1;3 -1 2;11 3 0];
>>B=[2 10 8]';
>>D=det(A)
>>[L,U]=lu(A)
>>X=U\\(L\\B)
显示结果如下：
D =
0
L =
0.3636 -0.5000 1.0000
0.2727 1.0000 0
1.0000 0 0
U =
11.0000 3.0000 0
0 -1.8182 2.0000
0 0 0.0000
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 2.018587e-017.
> In DMatlab\\pujun\\lx0720.m at line 4
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=[1 2 2 1;2 1 -2 -2;1 -1 -4 -3];
>>format rat %指定有理式格式输出
>>B=null(A,'r') %求解空间的有理基
运行后显示结果如下：
B =
2 5/3
-2 -4/3
1 0
0 1
或通过行最简行得到基：
>> B=rref(A)
B =
1.0000 0 -2.0000 -1.6667
0 1.0000 2.0000 1.3333
0 0 0 0
即可写出其基础解系（与上面结果一致）。
写出通解：
syms k1 k2
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 -2 3 -1;3 -1 5 -3;2 1 2 -2];
b=[1 2 3]';
B=[A b];
n=4;
R_A=rank(A)
R_B=rank(B)
format rat
if R_A==R_B&R_A==n %判断有唯一解
X=A\\b
elseif R_A==R_B&R_A<n %判断有无穷解
X=A\\b %求特解
C=null(A,'r') %求AX=0的基础解系
else X='equition no solve' %判断无解
end
运行后结果显示：
R_A =
2
R_B =
3
X =
equition no solve
说明 该方程组无解
例1-81 求解方程组的通解：
解法一：在Matlab编辑器中建立M文件如下：
A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
b=[1 4 0]';
B=[A b];
n=4;
R_A=rank(A)
R_B=rank(B)
format rat
if R_A==R_B&R_A==n
X=A\\b
elseif R_A==R_B&R_A<n
X=A\\b
C=null(A,'r')
else X='Equation has no solves'
end
运行后结果显示为：
R_A =
2
R_B =
2
Warning: Rank deficient, rank = 2 tol = 8.8373e-015.
> In DMatlab\\pujun\\lx0723.m at line 11
X =
0
0
-8/15
3/5
C =
3/2 -3/4
3/2 7/4
1 0
0 1
所以原方程组的通解为X=k1 +k2 +
解法二：用rref求解
A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
b=[1 4 0]';
B=[A b];
C=rref(B) %求增广矩阵的行最简形，可得最简同解方程组。
运行后结果显示为：
C =
1 0 -3/2 3/4 5/4
0 1 -3/2 -7/4 -1/4
0 0 0 0 0
对应齐次方程组的基础解系为： ， 非齐次方程组的特解为： 所以，原方程组的通解为：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。
>> load west0479
>> 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('iteration number') %x轴为迭代次数。
>> ylabel('relative residual') %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
>>load west0479;
>>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('iteration number')
>>ylabel('relative residual')
结果为
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*on 4*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)
gmres(A,b,restart,tol,maxit)
gmres(A,b,restart,tol,maxit,M)
gmres(A,b,restart,tol,maxit,M1,M2)
gmres(A,b,restart,tol,maxit,M1,M2,x0)
[x,flag] = gmres(A,b,…)
[x,flag,relres] = gmres(A,b,…)
[x,flag,relres,iter] = gmres(A,b,…)
[x,flag,relres,iter,resvec] = gmres(A,b,…)
除参数restart（缺省时为系数方阵A的阶数n）可以给出外，调用方式和返回的结果形式与命令bicg一样。
例1-86
>>load west0479;
>>A = west0479;
>>b = sum(A,2);
>>[x,flag] = gmres(A,b,5)
结果显示flag=1，表示在默认精度和默认迭代次数下不收敛。
若最后一行改为
>>[L1,U1] = luinc(A,1e-5);
>>[x1,flag1] = gmres(A,b,5,1e-6,5,L1,U1)
结果flag1=2，说明该方法失败，原因是U1的对角线上有0元素，构成坏条件的预处理因子。
若改为：
[L2,U2] = luinc(A,1e-6);
tol = 1e-15;
[x4,flag4,relres4,iter4,resvec4] = gmres(A,b,4,tol,5,L2,U2)
[x6,flag6,relres6,iter6,resvec6] = gmres(A,b,6,tol,3,L2,U2)
[x8,flag8,relres8,iter8,resvec8] = gmres(A,b,8,tol,3,L2,U2)
结果flag4=0，flag6=0，flag8=0表示参数restart分别取为4，6，8时收敛。
1.4.10 最小残差法解方程组
函数 minres
格式 x = minres(A,b)
minres(A,b,tol)
minres(A,b,tol,maxit)
minres(A,b,tol.maxit,M)
minres(A,b,tol,maxit,M1,M2)
minres(A,b,tol,maxit,M1,M2,x0)
[x,flag] = minres(A,b,…)
[x,flag,relres] = minres(A,b,…)
[x,flag,relres,iter] = minres(A,b,…)
[x,flag,relres,iter,resvec] = minres(A,b,…)
[x,flag,relres,iter,resvec,resveccg] = minres(A,b,…)
这里A为对称矩阵，这种方法是寻找最小残差来求x。调用方式和返回的结果形式与命令bicg一样。
例1-87
>>n = 100; on = ones(n,1);
>>A = spdiags([-2*on 4*on -2*on],-1:1,n,n);
>>b = sum(A,2);
>>tol = 1e-10;
>>maxit = 50;
>>M1 = spdiags(4*on,0,n,n);
>> [x,flag,relres,iter,resvec,resveccg] = minres(A,b,tol,maxit,M1,[],[])
结果显示：
flag =
0
relres =
4.6537e-014
iter =
49
resvec =（略）
resveccg =（略）
1.4.11 预处理共轭梯度方法
函数 pcg
格式 x = pcg(A,b)
pcg(A,b,tol)
pcg(A,b,tol,maxit)
pcg(A,b,tol,maxit,M)
pcg(A,b,tol,maxit,M1,M2)
pcg(A,b,tol,maxit,M1,M2,x0)
pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)
[x,flag] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)
[x,flag,relres] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)
[x,flag,relres,iter] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)
[x,flag,relres,iter,resvec] = pcg(A,b,tol,maxit,M1,M2,x0,p1,p2,…)
调用方式和返回的结果形式与命令bicg一样，这里A为对称正定矩阵。
1.4.12 准最小残差法解方程组
函数 qmr
格式 x = qmr(A,b)
qmr(A,b,tol)
qmr(A,b,tol,maxit)
qmr(A,b,tol,maxit,M)
qmr(A,b,tol,maxit,M1,M2)
qmr(A,b,tol,maxit,M1,M2,x0)
qmr(afun,b,tol,maxit,m1fun,m2fun,x0,p1,p2,…)
[x,flag] = qmr(A,b,…)
[x,flag,relres] = qmr(A,b,…)
[x,flag,relres,iter] = qmr(A,b,…)
[x,flag,relres,iter,resvec] = qmr(A,b,…)
调用方式和返回的结果形式与命令bicg一样，这里A为方阵。
例1-88
>>load west0479;
>>A = west0479;
>>b = sum(A,2);
>>[x,flag] = qmr(A,b)
结果显示flag=1，表示在缺省情况下不收敛
若改为
>>[L1,U1] = luinc(A,1e-5);
>>[x1,flag1] = qmr(A,b,1e-6,20,L1,U1)
结果显示 flag=2，表示是坏条件的预处理因子，不收敛。
若改为
>>[L2,U2] = luinc(A,1e-6);
>>[x2,flag2,relres2,iter2,resvec2] = qmr(A,b,1e-15,10,L2,U2)
>>semilogy(0:iter2,resvec2/norm(b),'-o') %每次迭代的相对残差图形
>>xlabel('iteration number')
>>ylabel('relative residual')
结果为
x=（全为1，略）
flag2 =
0 %表示收敛
relres2 = %表示相对残差
2.8715e-016
iter2 = %计算终止时的迭代次数
8
resvec2 = %每次迭代的残差
1.0e+005 *
7.0557
7.1773
3.4032
1.7220
0.0000
0.0000
0.0000
0.0000
0.0000
--------------------------------------------------------------------------------
特征值与二次型
1.5
工程技术中的一些问题，如振动问题和稳定性问题，常归结为求一个方阵的特征值和特征向量。
1.5.1 特征值与特征向量的求法
设A为n阶方阵，如果数“ ”和n维列向量x使得关系式 成立，则称 为方阵A的特征值，非零向量x称为A对应于特征值“ ”的特征向量。
详见1.3.5和1.3.6节：特征值分解问题。
例1-89 求矩阵 的特征值和特征向量
解：
>>A=[-2 1 1;0 2 0;-4 1 3];
>>[V,D]=eig(A)
结果显示：
V =
-0.7071 -0.2425 0.3015
0 0 0.9045
-0.7071 -0.9701 0.3015
D =
-1 0 0
0 2 0
0 0 2
即：特征值-1对应特征向量(-0.7071 0 -0.7071)T
特征值2对应特征向量(-0.2425 0 -0.9701)T和(-0.3015 0.9045 -0.3015)T
例1-90 求矩阵 的特征值和特征向量。
解：
>>A=[-1 1 0;-4 3 0;1 0 2];
>>[V,D]=eig(A)
结果显示为
V =
0 0.4082 -0.4082
0 0.8165 -0.8165
1.0000 -0.4082 0.4082
D =
2 0 0
0 1 0
0 0 1
说明 当特征值为1 (二重根)时，对应特征向量都是k (0.4082 0.8165 -0.4082)T，k为任意常数。
1.5.2 提高特征值的计算精度
函数 balance
格式 [T,B] = balance(A) %求相似变换矩阵T和平衡矩阵B，满足 。
B = balance(A) %求平衡矩阵B
1.5.3 复对角矩阵转化为实对角矩阵
函数 cdf2rdf
格式 [V,D] = cdf2rdf (v,d) %将复对角阵d变为实对角阵D，在对角线上，用2×2实数块代替共轭复数对。
例1-91
>> A=[1 2 3;0 4 5;0 -5 4];
>> [v,d]=eig(A)
v =
1.0000 -0.0191 - 0.4002i -0.0191 + 0.4002i
0 0 - 0.6479i 0 + 0.6479i
0 0.6479 0.6479
d =
1.0000 0 0
0 4.0000 + 5.0000i 0
0 0 4.0000 - 5.0000i
>> [V,D]=cdf2rdf(v,d)
V =
1.0000 -0.0191 -0.4002
0 0 -0.6479
0 0.6479 0
D =
1.0000 0 0
0 4.0000 5.0000
0 -5.0000 4.0000
1.5.4 正交基
命令 orth
格式 B=orth(A) %将矩阵A正交规范化，B的列与A的列具有相同的空间，B的列向量是正交向量，且满足：B'*B = eye(rank(A))。
例1-92 将矩阵 正交规范化。
解：
>>A=[4 0 0; 0 3 1; 0 1 3];
>>B=orth(A)
>>Q=B'*B
则显示结果为
P =
1.0000 0 0
0 0.7071 -0.7071
0 0.7071 0.7071
Q =
1.0000 0 0
0 1.0000 0
0 0 1.0000
1.5.5 二次型
例1-93 求一个正交变换X=PY，把二次型
化成标准形。
解：先写出二次型的实对称矩阵
在Matlab编辑器中建立M文件如下：
A=[0 1 1 -1;1 0 -1 1;1 -1 0 1;-1 1 1 0];
[P,D]=schur(A)
syms y1 y2 y3 y4
y=[y1;y2;y3;y4]；
X=vpa(P,2)*y %vpa表示可变精度计算，这里取2位精度
f=[y1 y2 y3 y4]*D*y
运行后结果显示如下：
P =
780/989 780/3691 1/2 -390/1351
780/3691 780/989 -1/2 390/1351
780/1351 -780/1351 -1/2 390/1351
0 0 1/2 1170/1351
D =
1 0 0 0
0 1 0 0
0 0 -3 0
0 0 0 1
X =
[ .79*y1+.21*y2+.50*y3-.29*y4]
[ .21*y1+.79*y2-.50*y3+.29*y4]
[ .56*y1-.56*y2-.50*y3+.29*y4]
[ .50*y3+.85*y4]
f =
y1^2+y2^2-3*y3^2+y4^2
即 f = |
|
/1 
IP卡
狗仔卡
发表于 2005-7-22 03:16:21
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡
楼主
发表于 2005-7-22 10:35:16