😺

2024/06/19に公開

# はじめに

(qiita に投稿していたものと完全に同じものです。）

# メイン

## サンプルコード

function get_square_root(n::Int, x0, RHS)
#nは摂動の最高次数 ε^(n-1)まで展開.
#RHSは知りたい平方根の値.
#x^2 = RHS+εにおいて、RHSはxが整数になるような数（平方数）でεが摂動項.x0^2=RHS.
#例えばx^2=5であれば、RHS=4, ε=1.
eps = RHS-x0^2
# シンボルの定義
ϵ, x = symbols("ϵ x", real=true)
syms = symbols("x:\$n", real=true)
eps_power = [ϵ^i for i in 0:n-1]

# x0 の設定とテイラー展開の形式で x を定義
x = sum([syms[i] * eps_power[i] for i in 1:n])
eq = x^2 - x0^2 - ϵ
f = expand(eq)

# ϵ の係数のリストを作成
coeff = Sym[[f.coeff(ϵ^i) for i in N(degree(f,gen=ϵ)):-1:1]; f(ϵ=>0)]

# 係数方程式を解くための変数を定義
solutions = collect(syms)

# for ループで係数方程式を解く
for i in 1:n
eq = coeff[end - i + 1]
for j in 1:i-1
eq = subs(eq, syms[j], solutions[j])
end
if i==1

solutions[i] = solve(eq, syms[i])[2]  #solve[2]が正の解、solve[1]が負の解
else
solutions[i] = solve(eq, syms[i])[1]
end
end

# 平方根の計算
eps_list = [eps^i for i in 0:length(solutions)]
x_sol = convert(Float64,sum([solutions[i] * eps_list[i] for i in 1:length(solutions)]))

# 結果の表示
return solutions, x_sol
end


from sympy import symbols, expand, solve

def get_square_root(n, x0, RHS):
# n is the highest order of perturbation, expand up to ε^(n-1).
# RHS is the value of the square root you want to know.
# In x^2 = RHS+ε, RHS is a number (square number) that makes x an integer, and ε is a perturbation term. x0^2=RHS.
# For example, if x^2=5, then RHS=4, ε=1.

eps = RHS - x0**2

# Define symbols
ϵ, x = symbols('ϵ x', real=True)
syms = symbols(f'x0:{n}', real=True)
eps_power = [ϵ**i for i in range(n)]

# Set x0 and define x in the form of Taylor expansion
x = sum([syms[i] * eps_power[i] for i in range(n)])
eq = x**2 - x0**2 - ϵ
f = expand(eq)

# Create a list of coefficients of ϵ
coeff = [f.coeff(ϵ, i) for i in reversed(range(n))] + [f.subs(ϵ, 0)]

# Define variables to solve the coefficient equation
solutions = [None] * n

# Solve the coefficient equation in a for loop
for i in range(n):
eq = coeff[n - i - 1]
for j in range(i):
eq = eq.subs(syms[j], solutions[j])
if i == 0:
solutions[i] = solve(eq, syms[i])[1]  # solve()[1] is the positive solution, solve()[0] is the negative solution
else:
solutions[i] = solve(eq, syms[i])[0]

# Calculate the square root
eps_list = [eps**i for i in range(len(solutions) + 1)]
x_sol = float(sum([solutions[i] * eps_list[i] for i in range(len(solutions))]))

# Display the result
return solutions, x_sol


\epsilon が負になる場合は試していないので、動作は保証しない。（欠陥コードだが、許してほしい。）
ただし、摂動論自体は \epsilon の符号について制限をかけてないので問題ないはず。

## 使い方

1. 関数の引数

• n\epsilonn-1 乗まで展開する
• x_0x_0^2<\text{(右辺)} である適当な数
• \text{RHS}：二次方程式の右辺
2. 出力
１つ目の出力（solutions）は係数方程式の解。出力方法はget_square_root(n,x0,RHS)[1]（pythonなら[1]\rightarrow[0]）
２つ目の出力（x_sol）は平方根の近似値。出力方法はget_square_root(n,x0,RHS)[2]（pythonなら[2]\rightarrow[1]）

### 例

\sqrt{5} を求めるときは、x^2 = 5 として正の解の近似値を求める。この場合、\text{RHS}=5 である。
コードは次のようになる：

get_square_root(20,2,5)


ただし、n は適当な整数で、x_0=2, x^2=4+\epsilon \ (\epsilon=1) とした。

# 興味深い(?)こと

get_square_root(5,2,5)[2]


この結果は、2.23602294921875 となった。

（量子力学の場合でも精度よくできるのか知っている人がいたら是非コメントなどで教えてください。）

ただし、\sqrt{2} を求めるときに、x_0=1, \ \epsilon=1とすると、精度が非常に悪くなる。

get_square_root(100,1,2)[2]


を実行すると、1.4143562059036923となる。
オーダーは同じなのに、\sqrt{5} に比べて精度が悪くなる理由は何なのだろうか。

pythonの動作確認をしているときに気づいたのだが、juliaより早く出力された。sympyとの相性だと勝手に思っているが、はたして......

# P.S.

また、問題設定は数式で書いた方が簡単なのだが、普通に面倒だったのでまた今度追記するかもしれない。