@ntsk2543posted at 2022-03-02

Pythonで漸化式を作成して配列に入れたい

Q&A

Closed

解決したいこと

漸化式の一般項が以下の様になっています。


f_n(x)=\frac{(x+r)^n+(x+r)^n}{2} \quad r=\sqrt{x^2-1}

この一般項を作成し$n$をfor文で回し、$x$についての関数にしたものをnumpy型の配列に保存したいです。
更に、保存したものを1つ取り出して、その関数に値を代入して得られた結果を別の配列に保存することが目標になります。

該当するソースコード

現在、最初の配列の代入で詰まっています。

import sympy as sp
import numpy as np

n=sp.Symbol("n",positive=True)
x=sp.Symbol("x",positive=True)


orbit=np.ones(100)
#一般項の作成(rは直接代入しています)
general_terms=((((x+sp.sqrt((x**2)-1))**n)+((x-sp.sqrt((x**2)-1))**n)))/2


for i in range(1,101):
    general_terms_2=general_terms.subs(n,i)
    orbit[i]=general_terms_2

発生している問題・エラー

TypeError: can't convert expression to float

また、$f_1(x)=x,f_2(x)=2x^2-1,f_3(x)=4x^3-3x$となるのですが、出力をみると

f1(x)=x
f2(x)=(x - sqrt(x**2 - 1))**2/2 + (x + sqrt(x**2 - 1))**2/2
f3(x)=(x - sqrt(x**2 - 1))**3/2 + (x + sqrt(x**2 - 1))**3/2

となっていて計算が簡単な式で表せていないことも問題になります。

自分で試したこと

lambdifyを用いてNumpyに変換してから代入をしてみました。

for i in range(1,101):
    general_terms_2=general_terms.subs(n,i)
    #n代入後はxの1変数のみになる。
    args = (x)
    func = sp.lambdify(args, general_terms_2, "numpy")
    orbit[i]=func

ただ、だめそうでした。

TypeError: float() argument must be a string or a number, not 'function'

0 likes

2Answer

@PondVillege posted at 2022-03-02

TypeError: can't convert expression to float

np.ones(100)で初期化したnumpy配列orbitでは，デフォルトの型がnumpy.float64で定義されており，orbit[i]=general_terms_2のようにして代入した場合はfloat型に変換するようにされているみたいです．

扱う型をsympy.core.symbol.Symbolに設定して初期化すれば動きます．

orbit=np.ones(100, dtype=sp.core.symbol.Symbol)

出力をみると計算が簡単な式で表せていない

計算を簡単な式で表すにはsympy.simplify()を使えばできます．
実際に$n$が1から10まで変化させて実行した結果を次に示します．

改訂案

import sympy as sp
import numpy as np

n=sp.Symbol("n",positive=True)
x=sp.Symbol("x",positive=True)


orbit=np.ones(100, dtype=sp.core.symbol.Symbol)
#一般項の作成(rは直接代入しています)
general_terms=((((x+sp.sqrt((x**2)-1))**n)+((x-sp.sqrt((x**2)-1))**n)))/2

for i in range(1,11):
    general_terms_2=general_terms.subs(n,i)
    print(sp.simplify(general_terms_2))
    orbit[i]=general_terms_2

出力結果

x
2*x**2 - 1
x*(4*x**2 - 3)
8*x**4 - 8*x**2 + 1
x*(16*x**4 - 20*x**2 + 5)
32*x**6 - 48*x**4 + 18*x**2 - 1
x*(64*x**6 - 112*x**4 + 56*x**2 - 7)
128*x**8 - 256*x**6 + 160*x**4 - 32*x**2 + 1
x*(256*x**8 - 576*x**6 + 432*x**4 - 120*x**2 + 9)
(x - sqrt(x**2 - 1))**10/2 + (x + sqrt(x**2 - 1))**10/2

もちろんシンプルな式になっているので，$f_3(x)$に関しては$f_3(x) = x(4x^2-3)$になります．

2Like

Comments

@ntsk2543
Questioner
これは,nが10以降の出力を見る限りympy.simplify()を用いても簡略化されない、ということはこれが一番簡単な出力であると考えられる、ということでよろしいのでしょうか？
もしくは、式が複雑(特にn乗の計算は大変)なので、表示しきれないということでしょうか？
@PondVillege
https://docs.sympy.org/latest/tutorial/simplification.html

ドキュメントを読むところ，heuristicに簡略化しているみたいで，「式が複雑だからこれ以上簡単にはできない」が正しそうです．2回sympy.simplify()を実行しても変化しませんでした．
@PondVillege
2回sympy.simplify()を実行しても変化しませんでしたが，2回目をsympy.factor()にしたら
(2*x**2 - 1)*(256*x**8 - 512*x**6 + 304*x**4 - 48*x**2 + 1)
になりました．
@ntsk2543
Questioner
非常に素早い回答、ありがとうございます。
ただ、私も失念しておりましたが質問の後半部分『更に、保存したものを1つ取り出して、その関数に値を代入して得られた結果を別の配列に保存することが目標になります。』について、ご教示いただけると非常に助かります。

@PondVillege posted at 2022-03-02

保存したものを1つ取り出して、その関数に値を代入して得られた結果を別の配列に保存する

orbitに保存されたもののうち第4項を取り出して，$1\le x \le 10$でx_sampleサンプル分の値を計算させました．

import sympy as sp
import numpy as np

n_sample, x_sample = 10, 100

n=sp.Symbol("n",positive=True)
x=sp.Symbol("x",positive=True)

orbit=np.ones(n_sample, dtype=sp.core.symbol.Symbol)
#一般項の作成(rは直接代入しています)
general_terms=((((x+sp.sqrt((x**2)-1))**n)+((x-sp.sqrt((x**2)-1))**n)))/2

for i in range(1, n_sample + 1):
    orbit[i]=sp.factor(sp.simplify(general_terms.subs(n,i)))

result = np.zeros((n_sample, x_sample))

result[3] = sp.lambdify(x, orbit[3], "numpy")(np.linspace(1, 10, x_sample))
print(result[3])

出てくる値はfloat型なので，np.zeros()の初期化の際には特に型指定をしていません．また，$x$の範囲を使いまわす場合は

x_values = np.linspace(1, 10, x_sample)
result[3] = sp.lambdify(x, orbit[3], "numpy")(x_values)

にしておいた方がよさそうです．複数取り出してグラフを書くなら

from matplotlib import pyplot as plt

x_values = np.linspace(1, 10, x_max)
for i in range(1, 5):
    result[i] = sp.lambdify(x, orbit[i], "numpy")(x_values)
    plt.plot(x_values, result[i], label=f"n={i}")

plt.xlabel("x")
plt.ylabel("y")
plt.yscale('log')
plt.grid()
plt.legend()
plt.show()

ですね．

1Like

Comments

@ntsk2543
Questioner
わざわざ丁寧に教えてくださり、ありがとうございました!

Are you sure you want to delete the question?

Pythonで漸化式を作成して配列に入れたい

解決したいこと

該当するソースコード

発生している問題・エラー

自分で試したこと

2Answer

Comments

Comments

Your answer might help someone💌