python コードのバグが分からない

下記のコードを書きました．
n=0,1,2,3,4,5までは上手くいくのですが
n=6以上に変えると計算がうまくいかないｔ，，，
とういうか計算時間が以上にかかるのですがなぜですか．
何かコードでおかしなことをしてますか？

たぶん，n=6以上にすると，
ret = solve(lp,x)で固まります．

#!/usr/bin/env python2.7
# coding: utf-8
# test.py

#変更する箇所は、ルジャンドルの欲しい次数（実際には欲しい零点の数)のn=だけ

import sympy
from sympy import Symbol, legendre, N, solve



list = []

x = Symbol('x')
n = 6
lp = legendre(n, x)
#print('多項式次数は', n)
#print(lp)  # 多項式を出力する

ret = solve(lp, x)  # lp = 0 を x について解く
#print('{}次ルジャンドル多項式の零点(解析解)は\n{}'.format(n, ret))

#解析解ではなく、数値として取り出す
vals = [N(r) for r in ret]
print('零点は')
#list.append(vals)
print(vals)

f1 = sympy.diff(lp,x)
#print('1階微分は')
#print(f1)

num = len(vals)

for i in range(num):
    f2 = f1.subs(x,vals[i])
    list.append(f2)
#print(list)

list2 = []

for i in range(num):
    w = 2/((1-vals[i]*vals[i])*list[i]*list[i])
    list2.append(w)
print('重みは')
print(list2)

firedfly

2019/01/20 04:14 編集

タイトルを見ただけで、なにをやろうとしてなにに困っているかわかるよう具体的に書きましょう。またコードについてもなにをするためのものなのか要約がほしいです。なぜだかわかりますか？回答者はまずタイトルを見ます。自分に経験があって回答できそうか、情報が多くて判断しやすいか。どんな試行錯誤を自分でしたかが書いておらず、長大なコードを丸投げして、ぜんぶ読んで実行して判断してくれという質問はもっとも負担の多いパターンです（そうしなければいけないこともありますが）。回答率を上げるため、そして自分の中で問題点を整理するため、なるたけ平易に最小限の説明・コードにしましょう。そのためには問題のスコープをできるだけ狭くすることです（問題になってるコードを特定して、それだけで実行できるようにしたものを載せる）。たとえば「sympyでルジャンドル多項式の次数を6以上にすると解くのに時間がかる」というタイトルと「Pythonコードのバグがわからない」ではだいぶ違うわけです（一般に計算時間がかかること自体はバグではないです）。コンピュータのリソースと同様に、回答者のリソースも有限です。なるたけ効率よく使えるようプログラム（質問）しましょう。

行動規範の内容に同意します

回答2件

firedflyさんが適切な回答をされました。
Tubasa1995さんのコメントで、
＞たぶん，n=6以上にすると，
＞ret = solve(lp,x)で固まります．
とされていますが、７分程で下記の結果がでました。

多項式次数は 6
231*x**6/16 - 315*x**4/16 + 105*x**2/16 - 5/16
6次ルジャンドル多項式の零点(解析解)は
[-(-395340976800*sqrt(3)*5**(1/6)*7**(1/3)*cos(atan(11*sqrt(6)/3)/3)**3/(2870175491568*5**(2/3)*7**(1/3)*sin(atan(11*sqrt(6)/3)/3)**4 + 5740350983136*5**(2/3)*7**(1/3)*sin(atan(11*sqrt(6)/3)/3)**2*cos(atan(11*sqrt(6)/3)/3)**2 + 2870175491568*5**(2/3)*7**(1/3)*cos(atan(11*sqrt(6)/3)/3)**4) - 158136390720*sqrt(3)*5**(2/3)*7**(1/3)*sin(2*atan(11*sqrt(6)/3)/3)*cos(atan(11*sqrt(6)/3)/3)**2/(2870175491568*5**(2/3)*7**(1/3)*sin(atan(11*sqrt(6)/3)/3)**4 + 5740350983136*5**(2/3)*7**(1/3)*sin(atan(11*sqrt(6)/3)/3)**2*cos(atan(11*sqrt(6)/3)/3)**2 + 2870175491568*5**(2/3)*7**(1/3)*cos(atan(11*sqrt(6)/3)/3)**4) - 395340976800*sqrt(3)*5
...省略
/2)/231]
数値解は
[-0.661209386466264 + 5.90432806037835e-103*I, 0.661209386466264 + 9.44692489660537e-102*I, -0.238619186083197 - 1.70461700631958e-102*I, 0.238619186083197 - 1.36369360505566e-101*I, -0.932469514203152 + 4.16328475581762e-103*I, 0.932469514203152 - 4.16328475581762e-103*I]
1階微分は
693*x**5/8 - 315*x**3/4 + 105*x/8
8.67837319736972 - 315*(0.661209386466264 + 9.44692489660537e-102*I)**3/4 + 1.23990889267945e-100*I + 693*(0.661209386466264 + 9.44692489660537e-102*I)**5/8
多項式次数は 2
3*x**2/2 - 1/2
2次ルジャンドル多項式の零点(解析解)は
[-sqrt(3)/3, sqrt(3)/3]
数値解は
[-0.577350269189626, 0.577350269189626]
1階微分は
3*x
1.73205080756888

待てば結果は出る事が確認されました。
固まる現象は出ませんでした。

投稿2019/01/20 10:07

oasis701

総合スコア148

ベストアンサー

こんにちは。

多項式を生成してsolverに任せているだけなので
コードにおかしいところは見当たりません。

内部的にどう解いているのかまでわからないですが
一般に多項式の次数を増やせば、解くのにかかる時間も増えます。
この単純な使い方でのsympyの限界、という認識でいいかと思います。

ちなみに、Mathmaticaのクラウド版ではさくっと解いてくれます。
231x**6/16 - 315x4/16 + 105*x2/16 - 5/16 = 0 - Wolfram|Alpha

下は参考用。
次数ごとに解くのにかかった時間を測定するプログラムです。
仰るとおり、n=6のときだけいきなり時間が急増していますね。

Python
1import time
2import sympy
3from sympy import Symbol, legendre, N, solve
4
5
6for n in range(1, 7):
7  x = Symbol('x')
8  lp = legendre(n, x)
9
10  start = time.time()
11  ret = solve(lp, x)
12  
13  elapse = round(time.time() - start, 3)
14  print("n:{}	time:{}	lp:{}".format(n, elapse, lp))

n:1	time:0.017	lp:x
n:2	time:0.026	lp:3*x**2/2 - 1/2
n:3	time:0.008	lp:5*x**3/2 - 3*x/2
n:4	time:0.071	lp:35*x**4/8 - 15*x**2/4 + 3/8
n:5	time:0.036	lp:63*x**5/8 - 35*x**3/4 + 15*x/8
n:6	time:248.679sec	lp:231*x**6/16 - 315*x**4/16 + 105*x**2/16 - 5/16

投稿2019/01/20 05:20

編集2019/01/20 05:27

firedfly

総合スコア1131