王笑权
本人为了学习过程的省心,在网上找“二次型化为标准型和规范型”的SAS程序,结果没有找到,便自己编写了一下。
基本思路是正交变换+初等变换结合。
暂时放在此地,省点本人的内存。
/*带有平方的二次型化解*/
proc iml;/*化为标准型的两步(1)原矩阵a求得特征向量e,(2)e`*a*e即可*/
a={17 -2 -2,-2 14 -4,-2 -4 14};/*二次型f=17x1**2+14x2**2-4x1x2-4x1x3-8x2x3例题来源 线性代数 P270/
田勇主编 机械工业出版社 2002.9 总策划 胡东华*/
call eigen(ei,e,a);/*特征函数*/
c=e`*a*e;/*特征向量与原矩阵积和情况就是正交变换要求得的标准二次型*/
c4=c//e;/*规范化第一步知e是由a的列出等变换而来*/
c3=c4[,1]#(1/sqrt(abs(c4[1])))||c4[,2]#(1/sqrt(abs(c4[2,2])))||c4[,3]#(1/sqrt(abs(c4[3,3])));/*用初等变换
化二次型为规范型第一步 接矩阵c4 对应列乘以原矩阵c的各列非0元素平方根倒数*/
c0=c3[4:6,1:3]`*a*c3[4:6,1:3];
print ei e c c4 c3 c0;/*输出特征值ei、特征向量e、对角阵c(标准二次型 根据化法不惟一)、转换为规范型的合并矩阵c4、
规范型矩阵c0(惟一,正负惯性指数位置安排可以随意) 其实就是惯性指数的位置和个数而已*/
/*不含平方项的二次型*/
proc iml;
a={0 1 1,1 0 -2,1 -2 0};/*例如2 化f=2x1x2+2x1x3-6x2x3成标准型并求所用的变换矩阵*/
ai1=a//i(ncol(a));/*初等列矩阵*/
ai1[1,]=ai1[1,]+ai1[2,];/*第2行加到第一行上*/
ai1[,1]=ai1[,1]+ai1[,2];/*第2列加到第一列上*/
ai1[2,]=ai1[2,]-ai1[1,]#ai1[2,1]/ai1[1];/*第21位化为0*/
ai1[3,]=ai1[3,]-ai1[1,]#ai1[3,1]/ai1[1];/*第31位化为0*/
ai1[,2]=ai1[,2]-ai1[,1]#ai1[1,2]/ai1[1];/*第12位化为0*/
ai1[,3]=ai1[,3]-ai1[,1]#ai1[1,3]/ai1[1];/*第22首位化为0*/
ai1[,3]=ai1[,3]-ai1[,2]#ai1[2,3]/ai1[2,2];/*第23位化为0*/
ai1[3,]=ai1[3,]-ai1[2,]#ai1[3,2]/ai1[2,2];/*第32位化为0 此时矩阵的上半部分就是化成的二次标准型(不惟一)*/
c1=ai1[4:6,];/*二次型转变为标准型的线性替换矩阵 就是ai1的下半部分*/
d1=c1`*a*c1;/*结果就是ai1的上半部分*/
ai2=ai1[,1]#(1/sqrt(abs(ai1[1])))||ai1[,2]#(1/sqrt(abs(ai1[2,2])))||ai1[,3]#(1/sqrt(abs(ai1[3,3])));
/*各列用上矩阵对角元素平方根倒数转化*/
c2=ai2[4:6,];/*二次型转变为规范型的线性替换矩阵 就是ai2的下半部分*/
d2=c2`*a*c2;/*结果规范型(惟一) 实质就是ai2的上半部分*/
print ai1 d1 d2;
结果:
SAS 系统 2008年04月11日 星期五 下午10时21分41秒 1
EI E C C4 C3
18 0.8944272 0.2981424 0.3333333 18 0 0 18 0 0 4.2426407 0 0
18 -0.447214 0.5962848 0.6666667 0 18 1.203E-15 0 18 1.203E-15 0 4.2426407 4.01E-16
9 0 -0.745356 0.6666667 0 1.145E-15 9 0 1.145E-15 9 0 2.699E-16 3
0.8944272 0.2981424 0.3333333 0.2108185 0.0702728 0.1111111
-0.447214 0.5962848 0.6666667 -0.105409 0.1405457 0.2222222
0 -0.745356 0.6666667 0 -0.175682 0.2222222
C0
1 0 0
0 1 1.4E-16
0 1.065E-16 1
SAS 系统 2008年04月11日 星期五 下午10时21分41秒 2
AI1 D1 D2
2 0 0 2 0 0 1 0 0
0 -0.5 0 0 -0.5 0 0 -1 0
0 0 4 0 0 4 0 0 1
1 -0.5 2
1 0.5 -1
0 0 1