如何提高 NumPy 中 FEM 计算的数字表示精度？

Question

为了我的硕士论文，我正在制作一个 FEM 程序，它需要很高的精度——我认为超过 30 位小数——但我不能这样做，因为 NumPy 的 float64 仅数据类型允许 16 个位置。

我试过 NumPy longdouble，但是当我尝试使用他们的 solve 方法时出现错误。

你们知道更简单的提高精度的方法吗？就像 MATLAB 甚至 Maple。提前致谢！

Answer 1

上述 FEM/HPC 从业者提出的几乎所有反对意见都是合理的，但已经接触过长期/低退化传热模拟（在大的、细粒度的时间尺度上）and/or 其他与数值处理相关的退化，主要来自 IEEE-754 表示的内在限制（在深度迭代等下）只需要超出标准数值处理的一步，
+ 我也尊重 OP 的需要/希望为她/他的进一步努力找到一些方向。

可行：_但成本...

在条件恶劣的数值分析中使用了以下方法，这很有用。

重新制定的成本函数和梯度函数以 pure_DEC 方式实现，接下来是 LSQ 最小化求解器。

问题不是本身的精度，而是一个聪明的问题-（重新）表述，以便利用高级 decimal-class 最小 [TIME] 域惩罚的内置操作。

>>> with decimal.localcontext() as locCTX:
       ...      for aPREC in range( 20, 31 ):
       ...          locCTX.prec = aPREC
       ...          ( pure_dec_LSQ_5DoF( locCTX, dec_fmin_x0_SEARCH_TRIM_TO_BE_PRECISE, decX, decY ), pure_dec_RESi( locCTX, dec_fmin_x0_SEARCH_TRIM_TO_BE_PRECISE, decX, decY ) )
       ...
       (Decimal('0.038471115298826195147'),           (Decimal('0.023589050081780503'),           Decimal('-0.082605913918299990'),           Decimal('0.150647690402532134'),           Decimal('-0.091630826566012630')))
       (Decimal('0.0384711152988261953165'),          (Decimal('0.0235890500817804889'),          Decimal('-0.0826059139182999933'),          Decimal('0.1506476904025321349'),          Decimal('-0.0916308265660126301')))
       (Decimal('0.03847111529882619531420'),         (Decimal('0.02358905008178048823'),         Decimal('-0.08260591391829999331'),         Decimal('0.15064769040253213501'),         Decimal('-0.09163082656601263007')))
       (Decimal('0.038471115298826195324048'),        (Decimal('0.023589050081780488368'),        Decimal('-0.082605913918299993309'),        Decimal('0.150647690402532135021'),        Decimal('-0.091630826566012630071')))
       (Decimal('0.0384711152988261953231489'),       (Decimal('0.0235890500817804883582'),       Decimal('-0.0826059139182999933087'),       Decimal('0.1506476904025321350199'),       Decimal('-0.0916308265660126300707')))
       (Decimal('0.03847111529882619532322276'),      (Decimal('0.02358905008178048835950'),      Decimal('-0.08260591391829999330863'),      Decimal('0.15064769040253213501998'),      Decimal('-0.09163082656601263007070')))
       (Decimal('0.038471115298826195323213788'),     (Decimal('0.023589050081780488359358'),     Decimal('-0.082605913918299993308625'),     Decimal('0.150647690402532135019974'),     Decimal('-0.091630826566012630070702')))
       (Decimal('0.0384711152988261953232136753'),    (Decimal('0.0235890500817804883593541'),    Decimal('-0.0826059139182999933086251'),    Decimal('0.1506476904025321350199740'),    Decimal('-0.0916308265660126300707023')))
       (Decimal('0.03847111529882619532321367314'),   (Decimal('0.02358905008178048835935336'),   Decimal('-0.08260591391829999330862505'),   Decimal('0.15064769040253213501997413'),   Decimal('-0.09163082656601263007070231')))
       (Decimal('0.038471115298826195323213665675'),  (Decimal('0.023589050081780488359353229'),  Decimal('-0.082605913918299993308625043'),  Decimal('0.150647690402532135019974132'),  Decimal('-0.091630826566012630070702306')))
       (Decimal('0.0384711152988261953232136649869'), (Decimal('0.0235890500817804883593532187'), Decimal('-0.0826059139182999933086250437'), Decimal('0.1506476904025321350199741307'), Decimal('-0.0916308265660126300707023064')))

从技术上讲，精度 -"expansion" 没有限制，但是，时间...:[=​​15=]

# [PERF] @ .prec ==    40               ~         4,000 [us] ~     4 [ms]
#        @ .prec == 10000 !!! --------- ~   991,875,234 [us] ~ 1,000 [ s]

这是为了获得预期的（主要是本身无限制）精度.[=15=而支付的成本]

import decimal
import numpy as np
import zmq

try:
    if  isinstance( decCTX, decimal.Context ):
        pass      # decCTX  already exists
except:
                    decCTX      = decimal.getcontext()
                    decCTX.prec = 60

decX = np.asarray( ( decimal.Decimal(  3.4 ), decimal.Decimal(  3.5 ), decimal.Decimal(  3.7 ), decimal.Decimal(   4.3 ), ) )
decY = np.asarray( ( decimal.Decimal( 65   ), decimal.Decimal( 85   ), decimal.Decimal( 97   ), decimal.Decimal( 100   ), ) )

...

dec_fmin_x0_SEARCH_ADAPTIVE = np.asarray( (                     decimal.Decimal( -101000000010553.05594055493064099456356276561617988943684402001075635                                                                                                                                                                                        ),
                                                                decimal.Decimal(               -8.660605201193546246                                                                                                                                                                                                                           ),
                                                                decimal.Decimal(                0.00021842459768549                                                                                                                                                                                                                            ),
                                                                decimal.Decimal(               99.9259163119085989057939988625810620201012857893012816197730189907743792931209843327426339987914746365315172977942868845721827684076717423116961495794648319380554868846324870276029626886129186998300662535940937605435069739237317269895772  ),
                                                                decimal.Decimal(                2.64971757369295002249999999827154484100152060917026952223212241653783649669777780217778380697777777796977777777969777777777969777777779697777777796977777805877778058777780587777777800577777780057777778005777777800577777774817774778285740 ),
                                                                )
                                           )