共形行列式

Question

你好，我有这个代码：

import numpy as np
from sympy import *
import sympy
 
x = Symbol('x')
x0 = Symbol('x0')
B = Symbol('B')
E = Symbol('E')
linewidth = Symbol('linewidth')

f = 1/pi*(linewidth/2)/((x - (E+2*B))**2 + (linewidth/2)**2)

def makeOneMeasurement(pCov,f,sigma,pvec,pmean):
    Iprior = np.linalg.inv(pCov)
    derivative_f = [f.diff(E),f.diff(B)]
    F = []
    for i in range(len(derivative_f)):
        for j in range(len(derivative_f)):
            g = derivative_f[i]*derivative_f[j]
            for k in range(len(pvec)):
                g = g.subs(pvec[k], pmean[k])
            F = F + [g]
    F = np.reshape(F,(2,2))/sigma**2
    F = sympy.Matrix(F)
    Iprior = sympy.Matrix(Iprior)
    
    return (Iprior + F).det()

如果我这样调用这个函数：

makeOneMeasurement([[1,0],[0,1]],f,0.1,["linewidth","E","B"],[0.1,990,0.11]).subs(x,990).evalf()

我得到输出1（正好是1）。但是，如果不返回return (Iprior + F).det()，我只返回return (Iprior + F)，我调用：

makeOneMeasurement([[1,0],[0,1]],f,0.1,["linewidth","E","B"],[0.1,990,0.11]).subs(x,990).evalf().det()

（基本上完全相同，但det在函数外部而不是内部调用），我得到：

3653.9564196053

我做错了什么？我希望能够将行列式作为x的函数返回（因此我希望使用第一个版本）。为什么我总是得到1而不是x的函数？谢谢大家!

Answer 1

更改代码以避免所有浮点数的使用，包括避免所有numpy例程的使用，如inv和reshape（使用sympy等价项）：

In [39]: import numpy as np
    ...: from sympy import *
    ...: import sympy
    ...: 
    ...: x = Symbol('x')
    ...: x0 = Symbol('x0')
    ...: B = Symbol('B')
    ...: E = Symbol('E')
    ...: linewidth = Symbol('linewidth')
    ...: 
    ...: f = 1/pi*(linewidth/2)/((x - (E+2*B))**2 + (linewidth/2)**2)
    ...: 
    ...: def makeOneMeasurement(pCov,f,sigma,pvec,pmean):
    ...:     Iprior = Matrix(pCov).inv()
    ...:     derivative_f = [f.diff(E),f.diff(B)]
    ...:     F = []
    ...:     for i in range(len(derivative_f)):
    ...:         for j in range(len(derivative_f)):
    ...:             g = derivative_f[i]*derivative_f[j]
    ...:             for k in range(len(pvec)):
    ...:                 g = g.subs(pvec[k], pmean[k])
    ...:             F = F + [g]
    ...:     F = Matrix(F).reshape(2,2)/sigma**2
    ...:     F = sympy.Matrix(F)
    ...:     Iprior = sympy.Matrix(Iprior)
    ...: 
    ...:     return (Iprior + F)
    ...: 

In [40]: M = makeOneMeasurement([[1,0],[0,1]],f,Rational('0.1'),["linewidth","E","B"],[Rational('0.1'
    ...: ),S(990),Rational('0.11')])

那么精确的矩阵和行列式是：

In [41]: M
Out[41]: 
⎡                   2                                          ⎤
⎢      ⎛      49511⎞              ⎛      49511⎞ ⎛      99022⎞  ⎥
⎢      ⎜2⋅x - ─────⎟              ⎜2⋅x - ─────⎟⋅⎜4⋅x - ─────⎟  ⎥
⎢      ⎝        25 ⎠              ⎝        25 ⎠ ⎝        25 ⎠  ⎥
⎢────────────────────────── + 1   ───────────────────────────  ⎥
⎢                         4                                 4  ⎥
⎢     ⎛           2      ⎞              ⎛           2      ⎞   ⎥
⎢   2 ⎜⎛    49511⎞     1 ⎟            2 ⎜⎛    49511⎞     1 ⎟   ⎥
⎢4⋅π ⋅⎜⎜x - ─────⎟  + ───⎟         4⋅π ⋅⎜⎜x - ─────⎟  + ───⎟   ⎥
⎢     ⎝⎝      50 ⎠    400⎠              ⎝⎝      50 ⎠    400⎠   ⎥
⎢                                                              ⎥
⎢                                                   2          ⎥
⎢ ⎛      49511⎞ ⎛      99022⎞          ⎛      99022⎞           ⎥
⎢ ⎜2⋅x - ─────⎟⋅⎜4⋅x - ─────⎟          ⎜4⋅x - ─────⎟           ⎥
⎢ ⎝        25 ⎠ ⎝        25 ⎠          ⎝        25 ⎠           ⎥
⎢ ───────────────────────────    ────────────────────────── + 1⎥
⎢                           4                             4    ⎥
⎢       ⎛           2      ⎞          ⎛           2      ⎞     ⎥
⎢     2 ⎜⎛    49511⎞     1 ⎟        2 ⎜⎛    49511⎞     1 ⎟     ⎥
⎢  4⋅π ⋅⎜⎜x - ─────⎟  + ───⎟     4⋅π ⋅⎜⎜x - ─────⎟  + ───⎟     ⎥
⎣       ⎝⎝      50 ⎠    400⎠          ⎝⎝      50 ⎠    400⎠     ⎦

In [42]: M.det()
Out[42]: 
                   2  8                         2  7                             2  6              
10000000000000000⋅π ⋅x  - 79217600000000000000⋅π ⋅x  + 274549981652000000000000⋅π ⋅x  - 54372976605
───────────────────────────────────────────────────────────────────────────────────────────────────
                                                        2  8                         2  7          
                                     10000000000000000⋅π ⋅x  - 79217600000000000000⋅π ⋅x  + 2745499

                  2  5                                   2  4                                      
8974880000000000⋅π ⋅x  + 673015111919049368767000000000⋅π ⋅x  - 533146420076341661747549392000000⋅π
───────────────────────────────────────────────────────────────────────────────────────────────────
                   2  6                                2  5                                   2  4 
81652000000000000⋅π ⋅x  - 543729766058974880000000000⋅π ⋅x  + 673015111919049368767000000000⋅π ⋅x  

2  3                      2                                         2  2                           
 ⋅x  + 50000000000000000⋅x  + 263966124524722600510236863978120000⋅π ⋅x  - 746812961137426061526108
───────────────────────────────────────────────────────────────────────────────────────────────────
                                     2  3                                         2  2             
- 533146420076341661747549392000000⋅π ⋅x  + 263966124524722600510236863978120000⋅π ⋅x  - 7468129611

                2                                                                                  
47700028830400⋅π ⋅x - 99022000000000000000⋅x + 49026782420000000000000 + 92438641532871794599827086
───────────────────────────────────────────────────────────────────────────────────────────────────
                              2                                               2                    
3742606152610847700028830400⋅π ⋅x + 9243864153287179459982708676348253060561⋅π                     

                2
76348253060561⋅π 
─────────────────
                 

正如你所看到的，里面有大量的数字。您可以使用精确替换获得行列式的公式，然后对其求值：

In [43]: M.det().subs(x, 990)
Out[43]: 
             2                   
67122964561⋅π  + 2420000000000000
─────────────────────────────────
                       2         
          67122964561⋅π          

In [44]: M.det().subs(x, 990).n()
Out[44]: 3653.95642136943

期望能够在不显著降低精度的情况下将低精度浮点替换为这样的表达式是不合理的，因此在最后调用evalf之前使用精确计算