mathematics计算矩阵特征值,奇异值
mathematica吧
全部回复
仅看楼主
level 6
2018年06月13日 13点06分 1
level 6
代码一直被百度吞楼,我想问的是,q-阶乘,用product,有问题吗
2018年06月13日 14点06分 8
level 6
Input["n=", n];
n = %
不删除,不吞楼
Input["q=", q];
q = %
B = IdentityMatrix[n + 1];
X = Array[x, n + 1];
X[[1]] = 1/(n + 2);
哈哈,今天天晴了
Do[X[[j]] = X[[1]] + (j - 1)/(n + 2), {j, 2, n + 1}]
f[r_, q_] := If[q == 1, r, (1 - q^r)/(1 - q)]
h[r_, q_] := Product[f[i, q], {i, 1, r}]
(*此处为新定义q-阶乘*)
2018年06月13日 14点06分 9
n=10,q=0.5
2018年06月13日 14点06分
level 6
接上面的,不善楼
g[r_, q_] := If[r == 0, 1, h[r, q]]
k[n_, i_, q_] :=
If[q == 1, Binomial[n, i], h[n, q]/(h[i, q] h[n - i, q])];
(*此处为定义的新组合数*)
z[q_, x_, n_, i_] :=
k[n, i, q]*x^i*Product[(1 - q^s*x), {s, 0, n - i - 1}]
Do[B[[1]][[j]] = z[q, X[[1]], n, j - 1], {j, 1, n + 1}]
哈哈哈,明天继续上课了
Do[B[[i]][[j]] = z[q, X[[i]], n, j - 1], {i, 2, n + 1}, {j, 1, n + 1}]
Eigenvalues[B] // N
这样球奇异值,特征值问题吗
NumberForm[SingularValueList[B, Tolerance -> 0], 30]
奇异值,如果只用sigular和特征值形式一样,,最后的两个求不出来,
2018年06月13日 14点06分 10
n=10,q=0.5
2018年06月13日 14点06分
level 6
t=
1.000000000000000e+00
4.702744772831446e-01
2.211116028093092e-01
1.192012063893431e-01
4.988050906080095e-02
1.503002970298066e-02
3.3
13343828318
723e-03
5.16
15809589478
32e-04
5.350997901013313e-05
3.293801296441902e-06
9.066819478722120e-08
这是matlab算法计算矩阵特征值,理论上与mathematica结果相对误差,在10的-15左右,现在与mathematics计算结果误差太大,算法没有问题,问题怀疑出在mathematica上
2018年06月13日 14点06分 11
level 6
2018年06月13日 14点06分 12
level 9
看你这从头造轮子还真是累。。。
mma内部实现了q函数,查看帮助文档guide/QFunctions
重新写了遍你的代码,事实上只需要这么几行就行
n=10;q=1/2;
z[q_, x_, n_, i_] := x^i QBinomial[n, i, q] QPochhammer[x, q, -i + n]
matrixB = Table[z[q, i/(n + 2), n, j - 1], {i, 1, n + 1}, {j, 1, n + 1}];
N[Eigenvalues[matrixB], 20]
N[SingularValueList[matrixB], 20]
B和其特征值、奇异值可以算出准确值
B的特征值数值近似后是
{1.0000000000000000000, 0.47027447728314387723, \
0.22111160280930921190, 0.11920120638934288330, \
0.049880509060800961035, 0.0
15030029702
980654701, \
0.0033
13343828318
7185466, 0.000516
15809589478
281592, \
0.000053509979010133102204, 3.2938012964419005842*10^-6,
9.0668194787221064636*10^-8}
B的奇异值数值近似后是
{1.62557529
13664580345
, 0.65496343042350221640, \
0.42614509301239471339, 0.16551650910525835661, \
0.048257582617247744837, 0.011984129075681537497, \
0.0023555777191550438126, 0.000351
18684456219
553458, \
0.00003735530
13752190019
66, 2.4867344053685638385*10^-6,
7.6095396928103957455*10^-8}
2018年06月14日 00点06分 14
谢谢你,我不会只好按照最笨的方法来,这里还想问一个问题,我发现了,q的输入,1/2,0.5造成了差别,我改成1/2后和你的结果就一样了,同样的若果把你的改成了0.5,结果就不是期待结果,这是为什呢?
2018年06月14日 07点06分
@爱吃大米😇 在mma中输入0.5会被视为机器精度数 在这里精度不够用导致了误差很大
2018年06月15日 01点06分
哦哦,也就是mathematica能用分数表示的数用分数是吧,另外我用DOT点乘了再计算特征值会有误差吗,DOT(matrixB,Transpose[matrixB])
2018年06月15日 02点06分
@爱吃大米😇 如果用精确数耗时可以接受 或者对精度要求十分高可以用 不然显式指定精确度也可以考虑下
2018年06月15日 02点06分
level 6
A = Dot[matrixB, Transpose[matrixB]];
N[Eigenvalues[A], 20]
N[SingularValueList[A], 20]
2018年06月15日 02点06分 15
你可以让mma显示出A。A每个元素都还是分数,也就是说没有精度损失。你还可以直接让mma计算Eigenvalues[A]和SingularValueList[A],你会得到使用Root函数表达的高次多项式的根,这也没有任何精度损失。N[Eigenvalues[A], 20]计算数值结果,20可以指定为任何你想要的精度(有效数字位数)。
2018年06月15日 02点06分
level 9
如图 假设你只要20位有效数字的精度
给定数字时只需要给比20多几位的精度就能与精确数计算相同的答案而时间节省了一半多(这里为了方便对计算过程计时改写成了一个模块)
关于数值计算精度的解释参见howto/ControlThePrecisionAndAccuracyOfNumericalResults
另外,只要计算中没有引入有限精度的数字且算法支持(一个不支持的例子的是FFT),mma默认进行准确数的运算
2018年06月15日 02点06分 16
看出来了,时间少了许多,谢谢您细心的指导,就是因为这个问题,我对matlab算法检查了许多遍,这样得到结果就是需要的结果
2018年06月15日 02点06分
1