python sympy グラフ描写

Question

最下部の画像の関数をpythonのsympyを使ってグラフ化しようと試みているのですが、計算時間が長いためかグラフが表示されず困っています。無限級数は数値結果を得やすくするために１０００までにしてあります。

どうにかして計算結果をグラフ化したいのですが、お教えいただけませんでしょうか？

初心者ですがよろしくお願いいたします。

m

Comments

[tiitoi]

コードはコピペできるように画像ではなく、 markdown 記法で記載していただけますか。

Answer 1

sympy は記号演算のライブラリなので、数式を記号で扱う上ではいいのですが、実際に値を計算しようとするとかなり遅いです。

グラフを描画する場合、具体的な数値を計算する必要があるので、numpy を使ったほうがいいでしょう。

計算する際は倍精度で行われますので、急激に値が増える階乗は!170あたりでオーバーフローします。
なので、1000項計算するのは数値計算の上では難しいでしょう。

numpy で計算する際のコツは for 文で1項ずつ計算して足していくのではなく、N 項分まとめて計算してから、和をとって級数にします。
matplotlib で描画する場合、描画範囲に離散的な点を例えば m 個作成して、その各点での関数値を計算する必要があるので、以下のようにします。

1. パラメータ t のときの各項 i の値を一度に計算する。形状 (m, N) の配列ができる。
2. 行方向に和をとる。形状 (m,) の配列ができる。

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import factorial, gamma

n = 100  # 計算する項数
l = 0.7  # ラムダ
t = np.linspace(0.01, 1, 100)  # 関数の描画範囲

k = np.arange(n + 1)

# 各項の符号
sign = np.where(k % 2 == 0, 1, -1)

# 各項の階乗の値
frac = factorial(k)

# 各項の sin 関数の値
sin = np.sin(np.pi * l * k)

# 各項のガンマ関数の値
g = gamma(l * k + 1)

# ガンマ関数の分母
t_ = t.reshape(-1, 1) ** (l * k + 1)

seq = sign / frac * sin * g / t_
print(seq.shape)  # (t, n) の配列

# 級数を求めるには、行方向に和を取ればよいので、
series = -1 / np.pi * seq.sum(axis=1)
print(series.shape)  # (100,)

fig, ax = plt.subplots()
ax.plot(t, series)

plt.show()

追記

1点質問なのですが、どうしても1000項以上計算したい場合は、どのような方法が考えられますでしょうか？アドバイスいただけますと幸いです。

1000項以上計算したい理由がグラフ描画であるならば、意味がないかもしれません。
例えば、頑張って有効数字20桁ぐらい計算したとしても、グラフ上ではその差は表現できないからです。
質問の式がなにかの関数のテイラー展開の結果なのだとしたら、最初の10項ぐらいでも有効数字が小数点以下何桁かはあると思いますが、グラフで描画する際の精度としては十分といえます。

その上でどうしても1000項ぐらい計算したい場合

一般的に扱われる倍精度 (Python の float 型) の範囲では、1000!は表現できる範囲を超えてしまうので、任意精度演算 のライブラリを使う必要があります。
任意精度演算ライブラリを使えば、メモリが許す限り、任意の桁数の値を扱うことができます。
ただし、欠点として、単精度/倍精度などハードウェアレベルで最適化されている演算ではないので、計算はかなり遅くなるでしょう。

Python で使える任意精度演算ライブラリとしては sympy についてくる mpmath があります。
float では170項ぐらいで階乗はオーバーフローしますが、2000ビットの多倍長演算で1000項まで計算でき、Wolfram Alpha の結果 と一致することを確認しました。

import sys

import mpmath
from scipy.special import factorial

mpmath.mp.dps = 2000  # 1000ビット

for i in range(1001):
    f = factorial(i)
    print(f"{i}! = {f:.2e}")

    if f > sys.float_info.max:
        break  # 171! で倍精度で表される範囲を超える。

f = mpmath.factorial(mpmath.mpf(1000))
print(f"1000! = {f}")

質問の式を mpmath で計算する例

import sys

import mpmath as mp

mpmath.mp.dps = 2000

s = 0
l = 0.3
n = 1000
t = 0.2

for k in range(n + 1):
    s += (
        (-1) ** k
        / mp.factorial(1000)
        * mp.sin(mpmath.pi * l * k)
        * mp.gamma(l * k + 1)
        / (t ** (l * k + 1))
    )

s *= -1 / mpmath.pi
print(s)
# -4.6852969663640794183052611442921398568086793342834616462334696324903867843231395630449988227880810131784478141659219945020203824927023753145340152275411511818373558920285872536714395228721948271826415724657244356679410653902566023614945014042952126178451891146305498687076853150690388109943827589984735451898842139557958621055625632539663091175942064908619685606545219293485810476695033253309084612857919257606784710540821141287519704017988233984166363710615914801845057844847498708240109737744243798479965501611309757531588120324018949625977650842536557835170122470528575372965177074318211049609908479423513835244125701516943697486654794175290551525218964214548834205190130209616768780024865572422562319780254088489151655831942096511712901359285285313156837232190787227348293663807776826088102389004523746146439727884614438552769772052583787278931259025862266267527347055478532365682291188279115875850275919795196649524369131042988665553942619715595168719062131399527261531325133906580769041387916239304444853504881787451700196905154750706687217882182421629624078044956529865543544068656579784986705519747832854686508939351378181527847396122106095160866115322720223585062715943851550247952235652437843682553103198759150039829159286341277372652777332702063251410719991844368057697345823111259489942513149539813031718828345330830157844803109246345868740928996257665817861367467687161910693184414807511579751981740167047404088833058340052830658962523612136648135596475773090234244305096492102145406041366097252467314623275357778005920563778105868762907482707819978490405043196844875157208281722161855859448349229516220667066123504490431317899702446327175215167205030224795654189487061289564293054352001359895654181422873855808898785005875909967936576766340843356222045846500038297018788320027808611619194276127030472392754256051166888548527686436583618833795960420300558879420743599986345799684146483437131704540648998092418292801233242752904307374935321733025575849004051292291211225721865580467562778e-1745

Comments

[kondasu]

お返事遅くなりすみません。丁寧なご回答ありがとうございます。
1点質問なのですが、どうしても1000項以上計算したい場合は、どのような方法が考えられますでしょうか？アドバイスいただけますと幸いです。

[tiitoi]

追記しました。
どうしても1000項計算したい場合、任意精度演算ライブラリをお使いください。
例えば、Python では sympy についてくる mpmath があります。
numpy や math モジュールの関数は使えないので、計算は全部 mpmath の関数で行う必要があります。
https://docs.sympy.org/0.7.1/modules/mpmath/functions/gamma.html

[tiitoi]

第1000項まで計算する例も記載しました。
第1000項まで計算するのに、1分ぐらいかかったので、t を [0.01, 1] の範囲で100個だとすると、約2時間近く計算にかかります。