摂動と平方根の値

はじめに

(qiita に投稿していたものと完全に同じものです。）
物理をやっている、あるいはやっていた人なら量子力学で摂動を見たことがあるだろう。この記事では、摂動を理解するための補助輪として、平方根の値を求めるコードを書いてみた。

メイン

サンプルコード

下のコードは Julia のもの。

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

以下、python 版。Julia のコードを GitHub Copilot で書き換えたもの。

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$：$\epsilon$ の $n-1$ 乗まで展開する
   • $x_0$：$x_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)$ とした。

興味深い(?)こと

上の例において、$n=5$ にしてみたのが次のコードである。

get_square_root(5,2,5)[2]

この結果は、2.23602294921875 となった。
量子力学の摂動論では $\epsilon$ が微小量と書かれているが、平方根の値を求める文脈においてはオーダーが一緒でも十分精度が良い。
（量子力学の場合でも精度よくできるのか知っている人がいたら是非コメントなどで教えてください。）

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

get_square_root(100,1,2)[2]

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

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

P.S.

摂動論については量子力学でなくても様々なところで使われている。ここでは、それを2次方程式に使っただけである。
また、問題設定は数式で書いた方が簡単なのだが、普通に面倒だったのでまた今度追記するかもしれない。