How to increase a number-representation precision for a FEM comptuting in NumPy?

Question

For my master's dissertation, I am making a FEM program that requires a lot of precision — more than 30 decimal places, I think — but I can't do that as NumPy's float64 data-type only allows 16 places.

I've tried NumPy longdouble, but it gives me an error when I try to use their solve method.

Do you guys know a simpler way to increase the precision? Just like MATLAB or even Maple. Thanks in advance!

Have you tried `decimal` module? Or is `numpy` a requirement? — jpp, Apr 24 '18 at 09:08
numpy is a requirement, as i do a lot of operations with matrices and arrays. But can i use let's say, for example, to increase the precision of a dot product? — tomas-silveira, Apr 24 '18 at 09:10
There is no support in numpy for more than 80bit float [longdouble](https://docs.scipy.org/doc/numpy-1.14.0/user/basics.types.html). What error do you get ? — Eolmar, Apr 24 '18 at 09:17
I'm getting "TypeError: array type float64 is unsupported in linalg" — tomas-silveira, Apr 24 '18 at 09:19
can you give a [mcve](https://stackoverflow.com/help/mcve)? `numpy.linalg` work with float64 on my machine — Eolmar, Apr 24 '18 at 09:25
Are you sure that precision is the problem? FEM is an approximation method,. I would review the theory before trying to work with bigger floats. — hpaulj, Apr 24 '18 at 09:33
It is working indeed on mine too Eolmar! Just a typo in my code. The code is too big to post ( > 1k lines code). I've reviewed the theory for the last 2 days with my professor, and all the functions are working properly (for a small number of elements), but when i try the program with more elements, it gives me bad results, which made us believe the problem might have something to do with precision. — tomas-silveira, Apr 24 '18 at 10:06
I couldn't agree more with hpaulj: a precision of 30 digits seems to be much higher than anything FEM or any other method could provide. And, as Eolmar suggested, creating a mcve could be in itself a great tool to analyse what's wrong with your code and find a solution, as well as some unit tests. Good luck with your work! — Gianluca Micchi, Apr 24 '18 at 12:23
It might help if you gave the usual FEM statistics - number of elements , nodes, nodes per element, basis functions, dof, size of the stiffness matrix etc. Also are you using just numpy, or scipy sparse? What solver? — hpaulj, Apr 24 '18 at 15:55
`numpy.linalg` and `scipy.sparse.linalg` solvers use fast compiled libraries where possible. As a result we are limited to the float types that those libraries implement. — hpaulj, Apr 24 '18 at 20:10

score 0 · Accepted Answer · edited May 11 '18 at 13:02

Almost all objections presented from FEM / HPC practitioners above are legitimate, yet having been exposed to cases, where long-term / low-degradation heat-transfer simulations ( on a large, fine-grained time-scale ) and/or other numerical-processing related degradations, principally coming from intrinsic limits of IEEE-754 representations ( under deep re-iterations et al ) simply required a step out of the standard numerical processing,
+ I also respect the OP's need / wish to see some direction for her / his further efforts.

Doable : _{yet at cost ...}

Used the following approach in badly-conditioned numerical analyses, where this worked as charm.

Reformulated cost-function and gradient-function were implemented in pure_DEC-fashion, next the LSQ-minimisation solver.

The problem is not the precision per-se, but a smart problem-(re)formulation, so as to harness the advanced decimal-class built-in operations at minimum [TIME]-domain penalty.

>>> 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')))

Technically, the precision-"expansion" goes without limits, yet, the time...:

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

This is the cost to pay so as to receive the expected ( principally unlimited per-se ) precision.

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 ),
                                                                )
                                           )

Thanks for your solution! altough i already found a workaround for my problem, this might be useful for the future when i need more decimal cases :) — tomas-silveira, Apr 25 '18 at 17:38

How to increase a number-representation precision for a FEM comptuting in NumPy?

1 Answers1

Doable : yet at cost ...

Doable : _{yet at cost ...}