注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Mr.Right

不顾一切的去想,于是我们有了梦想。脚踏实地的去做,于是梦想成了现实。

 
 
 

日志

 
 
关于我

人生一年又一年,只要每年都有所积累,有所成长,都有那么一次自己认为满意的花开时刻就好。即使一时不顺,也要敞开胸怀。生命的荣枯并不是简单的重复,一时的得失不是成败的尺度。花开不是荣耀,而是一个美丽的结束,花谢也不是耻辱,而是一个低调的开始。

网易考拉推荐

阿英讲二次型标准化 --- maple实现  

2013-08-18 15:11:55|  分类: 编程 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

即以此功德,庄严佛净土。上报四重恩,下济三途苦。惟愿见闻者,悉发菩提心。在世富贵全,往生极乐国。
# change quadratic form into standard diagonalized form 
restart;
with(LinearAlgebra):
with(VectorCalculus):
# example 1: quadratic form in 2 variables
Q1 := -6*x^2 + 4*x*y - 3*y^2;
H1 := Hessian(Q1, [x, y]);
D1 := DiagonalMatrix(Eigenvalues(H1)); 
P1 := IsSimilar(H1, D1, output ='C'); 
Equal(P1 . H1 . MatrixInverse(P1) , D1) ; # evalm(P1 &* H1 &* MatrixInverse(P1)) = D1
V1 := Eigenvectors(H1); 
V1 := op(2, [V1]);
# We would like to use the inverse of  P1 to define our change of variables
# Q, R := QRDecomposition(V1, fullspan); # 线性变换在不同基下的两个矩阵是相似的
# 右特征矢量V1的每一列实际上是“老基{x,y}”在“新基{u, v}”下的坐标!MatrixInverse(P1)和V1的效果是相同的
Q, R := QRDecomposition(MatrixInverse(P1), fullspan);
T := evalm(Q &* Matrix(2, 1, [u, v]));
out1 := subs(x = T[1, 1], y = T[2, 1], Q1);
printf("the coefficients are twice the eigenvalues of H1\n");
"standard quadratic form" = expand(out1);
阿英讲二次型标准化 --- maple实现 - 阿英 - Mr.Right
 
# example 2 : quadratic form in 3 variables
Q2 := 3*x^2 + 12*x*y +24*y*z - 6*x*z + 6*y^2 + 6*z^2;
H2 := Hessian(Q2, [x, y, z]);
D2 := DiagonalMatrix(Eigenvalues(H2)); 
P2 := IsSimilar(H2, D2, output ='C'); 
Equal(P2 . H2 . MatrixInverse(P2) , D2) ; # checking if equal
MatrixNorm(P2.H2.MatrixInverse(P2)-D2, Frobenius);
V2 := Eigenvectors(H2); 
V2 := op(2, [V2]);
# We would like to use the inverse of  P2 to define our change of variables
Q, R := QRDecomposition(MatrixInverse(P2), fullspan);
T := evalm(Q &* Matrix(3, 1, [u, v, w]));
out2 := subs(x = T[1, 1], y = T[2, 1], z = T[3, 1], Q2);
"standard quadratic form" = expand(out2);
# It would be nice if we could collect the terms as well
printf("the coefficients are the eigenvalues of H2\n");
result := collect(expand(out2), [u, v, w]):
"out2 after collect" = result;
printf("The coefficients do not look like twice the original eigenvalues, so let's check them\n");
coeff1 := subs( u = 1, v = 0, w = 0, result);
is(2*coeff1 = D2[1, 1]);
coeff2 := subs( u = 0, v = 1, w = 0, result);
is(2*coeff2 = D2[2, 2]);
coeff3 := subs( u = 0, v = 0, w = 1, result);
is(2*coeff3 = D2[3, 3]);

---------------------------------------
# 配方法
restart;
with(LinearAlgebra):
with(VectorCalculus):
# constraint 2
x := Matrix(2, 1, [u, v]);
A := Matrix([[1/4, 0], [0, 1]]);
b := Matrix(2, 1, [1/2, 0]);
y := Transpose(x) . A . x - Transpose(b) . x + 1/4;
Student[Precalculus][CompleteSquare](y[1, 1], u, v);


# constraint 3
x := Matrix(2, 1, [u, v]);
A := Matrix([[5/8, 3/8], [3/8, 5/8]]);
b := Matrix(2, 1, [11/2, 13/2]);
y := Transpose(x) . A . x - Transpose(b) . x + 35/2 + 1;
Student[Precalculus][CompleteSquare](y[1, 1], u, v);

  评论这张
 
阅读(828)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016