最下部の画像の関数をpythonのsympyを使ってグラフ化しようと試みているのですが、計算時間が長いためかグラフが表示されず困っています。無限級数は数値結果を得やすくするために1000までにしてあります。
どうにかして計算結果をグラフ化したいのですが、お教えいただけませんでしょうか?
初心者ですがよろしくお願いいたします。
気になる質問をクリップする
クリップした質問は、後からいつでもMYページで確認できます。
またクリップした質問に回答があった際、通知やメールを受け取ることができます。
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。

回答1件
0
ベストアンサー
sympy は記号演算のライブラリなので、数式を記号で扱う上ではいいのですが、実際に値を計算しようとするとかなり遅いです。
グラフを描画する場合、具体的な数値を計算する必要があるので、numpy を使ったほうがいいでしょう。
計算する際は倍精度で行われますので、急激に値が増える階乗は!170あたりでオーバーフローします。
なので、1000項計算するのは数値計算の上では難しいでしょう。
numpy で計算する際のコツは for 文で1項ずつ計算して足していくのではなく、N 項分まとめて計算してから、和をとって級数にします。
matplotlib で描画する場合、描画範囲に離散的な点を例えば m 個作成して、その各点での関数値を計算する必要があるので、以下のようにします。
- パラメータ t のときの各項 i の値を一度に計算する。形状 (m, N) の配列ができる。
- 行方向に和をとる。形状 (m,) の配列ができる。
python
1import numpy as np 2import matplotlib.pyplot as plt 3from scipy.special import factorial, gamma 4 5n = 100 # 計算する項数 6l = 0.7 # ラムダ 7t = np.linspace(0.01, 1, 100) # 関数の描画範囲 8 9k = np.arange(n + 1) 10 11# 各項の符号 12sign = np.where(k % 2 == 0, 1, -1) 13 14# 各項の階乗の値 15frac = factorial(k) 16 17# 各項の sin 関数の値 18sin = np.sin(np.pi * l * k) 19 20# 各項のガンマ関数の値 21g = gamma(l * k + 1) 22 23# ガンマ関数の分母 24t_ = t.reshape(-1, 1) ** (l * k + 1) 25 26seq = sign / frac * sin * g / t_ 27print(seq.shape) # (t, n) の配列 28 29# 級数を求めるには、行方向に和を取ればよいので、 30series = -1 / np.pi * seq.sum(axis=1) 31print(series.shape) # (100,) 32 33fig, ax = plt.subplots() 34ax.plot(t, series) 35 36plt.show()
追記
1点質問なのですが、どうしても1000項以上計算したい場合は、どのような方法が考えられますでしょうか?アドバイスいただけますと幸いです。
1000項以上計算したい理由がグラフ描画であるならば、意味がないかもしれません。
例えば、頑張って有効数字20桁ぐらい計算したとしても、グラフ上ではその差は表現できないからです。
質問の式がなにかの関数のテイラー展開の結果なのだとしたら、最初の10項ぐらいでも有効数字が小数点以下何桁かはあると思いますが、グラフで描画する際の精度としては十分といえます。
その上でどうしても1000項ぐらい計算したい場合
一般的に扱われる倍精度 (Python の float 型) の範囲では、1000!は表現できる範囲を超えてしまうので、任意精度演算 のライブラリを使う必要があります。
任意精度演算ライブラリを使えば、メモリが許す限り、任意の桁数の値を扱うことができます。
ただし、欠点として、単精度/倍精度などハードウェアレベルで最適化されている演算ではないので、計算はかなり遅くなるでしょう。
Python で使える任意精度演算ライブラリとしては sympy についてくる mpmath があります。
float では170項ぐらいで階乗はオーバーフローしますが、2000ビットの多倍長演算で1000項まで計算でき、Wolfram Alpha の結果 と一致することを確認しました。
python
1import sys 2 3import mpmath 4from scipy.special import factorial 5 6mpmath.mp.dps = 2000 # 1000ビット 7 8for i in range(1001): 9 f = factorial(i) 10 print(f"{i}! = {f:.2e}") 11 12 if f > sys.float_info.max: 13 break # 171! で倍精度で表される範囲を超える。 14 15f = mpmath.factorial(mpmath.mpf(1000)) 16print(f"1000! = {f}")
質問の式を mpmath で計算する例
python
1import sys 2 3import mpmath as mp 4 5mpmath.mp.dps = 2000 6 7s = 0 8l = 0.3 9n = 1000 10t = 0.2 11 12for k in range(n + 1): 13 s += ( 14 (-1) ** k 15 / mp.factorial(1000) 16 * mp.sin(mpmath.pi * l * k) 17 * mp.gamma(l * k + 1) 18 / (t ** (l * k + 1)) 19 ) 20 21s *= -1 / mpmath.pi 22print(s) 23# -4.6852969663640794183052611442921398568086793342834616462334696324903867843231395630449988227880810131784478141659219945020203824927023753145340152275411511818373558920285872536714395228721948271826415724657244356679410653902566023614945014042952126178451891146305498687076853150690388109943827589984735451898842139557958621055625632539663091175942064908619685606545219293485810476695033253309084612857919257606784710540821141287519704017988233984166363710615914801845057844847498708240109737744243798479965501611309757531588120324018949625977650842536557835170122470528575372965177074318211049609908479423513835244125701516943697486654794175290551525218964214548834205190130209616768780024865572422562319780254088489151655831942096511712901359285285313156837232190787227348293663807776826088102389004523746146439727884614438552769772052583787278931259025862266267527347055478532365682291188279115875850275919795196649524369131042988665553942619715595168719062131399527261531325133906580769041387916239304444853504881787451700196905154750706687217882182421629624078044956529865543544068656579784986705519747832854686508939351378181527847396122106095160866115322720223585062715943851550247952235652437843682553103198759150039829159286341277372652777332702063251410719991844368057697345823111259489942513149539813031718828345330830157844803109246345868740928996257665817861367467687161910693184414807511579751981740167047404088833058340052830658962523612136648135596475773090234244305096492102145406041366097252467314623275357778005920563778105868762907482707819978490405043196844875157208281722161855859448349229516220667066123504490431317899702446327175215167205030224795654189487061289564293054352001359895654181422873855808898785005875909967936576766340843356222045846500038297018788320027808611619194276127030472392754256051166888548527686436583618833795960420300558879420743599986345799684146483437131704540648998092418292801233242752904307374935321733025575849004051292291211225721865580467562778e-1745
投稿2019/08/16 11:34
編集2019/08/19 04:27総合スコア21960
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。

あなたの回答
tips
太字
斜体
打ち消し線
見出し
引用テキストの挿入
コードの挿入
リンクの挿入
リストの挿入
番号リストの挿入
表の挿入
水平線の挿入
プレビュー
質問の解決につながる回答をしましょう。 サンプルコードなど、より具体的な説明があると質問者の理解の助けになります。 また、読む側のことを考えた、分かりやすい文章を心がけましょう。