パイソンで計算の精度を上げる方法を教えてください。

Question

前提・実現したいこと

パイソンで計算の精度を上げる方法を教えてください。
今回、流体の数値解析でパイソンを使って計算しています。
初期の値は対象的な数値なのですが
[0.6366197723675814, -0.6366197723675814]
計算後には
[0.4379113806521834, -0.4467526621560704]
値としては近いのですが、対象的ではない値になってしまいました。
（対象的な値になるように計算しています。）

計算課程でnumpy.sin, cos, numpy.dotで2x６行列などを途中で計算しており、
検索して見たところパイソンでfloat計算すると計算誤差が多く出てしまうという記事をいくつも見つけました。
一番解決方法で多かった、Decimal と Fraction を使った計算も試して見たのですが、sin, cos の計算など、を含めるとなかなか精度が上げられず、期待した答えを出すことができませんでした。

どなたか、パイソンで計算精度を上げる方法など教えていただけないでしょうか。

主に、

　　　　　　　a0L=a0[i-1] if (i > 0) else a0[-1]
            a1L=a1[i-1] if (i > 0) else a1[-1]

            a0R=a0[i+1] if (i < Ne-1) else a0[0]
            a1R=a1[i+1] if (i < Ne-1) else a1[0]

            #Put current potision
            a0C=a0[i]
            a1C=a1[i]

            a_n = np.array([ [a0C],[a1C] ])
            a_tilda = np.array( [[a0L],[a1L],[a0C],[a1C],[a0R],[a1R]] )

         
            a_1 = a_n + (1/2)*D*dt/(h**2) * np.dot(A,a_tilda)
            a1_tilda = np.array( [[a0L],[a1L],a_1[0],a_1[1],[a0R],[a1R]] )
            a_n_1 = a_n + (dt/2)*(D/h**2)*np.dot(A,a_tilda)+ (dt/2)*(D/h**2)*np.dot(A,a1_tilda)

の計算課程で、計算精度が落ちているように感じております。a1, a0 は全て　float 型です。

a0[i] = ( np.exp( -D*((2*math.pi/L)**2)*t ) ) * ( (L/(h*math.pi)) *np.sin((h*math.pi)/L) ) * np.sin( (2*math.pi/L)*x[i] )

a1[i] = ( 3*np.exp( -D*((2*math.pi/L)**2)*t ) ) * (L/(h*math.pi))* ( ( L/(h*math.pi))*np.sin((h*math.pi)/L) - np.cos((h*math.pi)/L))*(np.cos( (2*math.pi/L)*x[i] ))

Decimal型にできるところは変更しても、求めたい数値が出ずに困っております。
どなたか、numpyを使っても、複雑な計算問題でも使える、精度を上げる方法を教えていただけないでしょうか。
よろしくお願いいたします。

Comments

[ikapy]

「Pythonでは decimal というパッケージを使えば任意精度の10進演算ができます」という記事をみつけました。
これを利用しては・・・

[ikapy]

https://oku.edu.mie-u.ac.jp/~okumura/python/decimal.html
です。

[alpha_high]

説明不足で申し訳ありません。
numpyを多数使用しており、できればnumpyを使用しつつDecimalを計算したいと思っております。
このような場合は、numpyを諦めて、行列計算やsin,cosなども.itemか何かで数字を呼び出してDecimalでコツコツ計算する、という方法しかないのでしょうか。かなり長い行数のコードの為、なるべく形は残しつつ精度を上げたいのですが、pythonにはDecimalしかないでしょうか。

[ikapy]

>pythonにはDecimalしかないでしょうか。

わかりませんです。
整数は多倍精度らしいけど、floatが多倍精度だと聞いた事は（今のところ）ありません。僕はご利益（制度）が欲しいなら、お布施（努力）も仕方ないかと思うけど。欲しいものがあることを！

[ikapy]

制度は精度ね

[ikapy]

わかりませんです。といいながら恥をかきますが、
「多倍長評価でnumpy.ndarrayを使う」というのを見つけました。
https://qiita.com/hanaata/items/7e466acb8897e6eac8d6
です。
ご所望のものかどうかは質問者さんが判断してください。

[jbpb0]

> a1, a0 は全て　float 型です。

当方のMacで確認したら、「float128」にすると、ちょっと計算結果が変わりました

import numpy as np
np.set_printoptions(precision=20)

a = np.array([1, 2, 3], dtype=np.float)
print(a.dtype)
aa = np.sin(a)
print(aa.dtype)
print(aa)

b = np.array([1, 2, 3], dtype=np.float128)
print(b.dtype)
bb = np.sin(b)
print(bb.dtype)
print(bb)

print(aa - bb)

差はほんのちょっとなので、質問者さんのケースで効果あるかは分かりませんが

[alpha_high]

ikapy さん
リンクをありがとうございます。Decimalでもnumpay.ndarrayが使えるのですね。。
sincosもループでなんとかすればDecimalで出せるという記事も見つけたので、検討してみます。本当にありがとうございます。

Jbpb0 さん
私のマックでもちょっと結果が変わってました。float128で書き直してみようと思います。ありがとうございます。

Answer 1

• パイソンでfloat計算すると計算誤差が多く出てしまうという記事をいくつも見つけました。

ウソではありませんが、誤解を招きやすい表現です。
現在のコンピュータでfloat計算すると、FORTRANであれ、C/C++であれ、Pythonであれ、その他どんなプログラミング言語を使っても、計算誤差が多く出てしまうことがあります。
これは全てのディジタルコンピュータハードウェアの性質ですので、プログラミング言語を変更しても同じように計算誤差が出ます。

Decimalを使うのは、コンピュータの内部で十進表現を使うことで、二進表現と人間普通使っている十進表現の差を埋めるためであって数値計算の誤差を減らすためには役に立ちません。

計算誤差を減らすためには、コンピュータのアーキテクチャ(特に浮動小数点レジスタの格納形式、演算方法)を理解し、誤差が起こる原因を理解し、誤差の原因に応じて演算順序の変更、共通項の括り出し、恒等式の変換などのテクニックを使ってソースコードを変更していきます。

少し検索したところでは、秋田大学 プログラミングテクニック（計算誤差）が比較的良くまとまっていましたので、これを熟読することをお勧めします。ここに書かれているのはC言語ですが、FORTRANでもPythonでもやることは同じです。

Comments

[ikadzuchi]

decimalは任意精度ですので数値計算の誤差を減らせます。