자동 미분 & 비선형 역학계 수치 계산

비선형 역학 분야에서 고정점의 안정성을 논의할 때 선형 안정성 분석을 한다.예를 들어 Brusselator 방정식를 예로 들면 그 역학 시스템은 다음과 같은 공식으로 표시한다($a, b$는 정상수).

\begin{equation}
\begin{aligned}
&\dot{x} = f(x,y) = a-(b+1)x+x^2y\\
&\dot{y} = g(x,y) = bx-x^2y
\end{aligned},\qquad x,y\ge 0
\end{equation}

이 방정식의 고정점이 $(x, y)=(x^*, y^*)$이면 선형 안정성 분석을 통해 고정점 $(x^*, y^*)$의 안정성을 계산하여 다음 아크비안 행렬 공식의 대각과 토론할 수 있다.

\begin{equation}
\frac{\partial (f,g)}{\partial (x,y)}=
\left(\begin{aligned}
&\frac{\partial f(x,y)}{\partial x}&\frac{\partial f(x,y)}{\partial y}\\
&\frac{\partial g(x,y)}{\partial x}&\frac{\partial g(x,y)}{\partial y}\\
\end{aligned}\right)_{\begin{aligned}x=x^*\\y=y^*\end{aligned}}
\end{equation}

아래의 그림은 고정점 중의 대각과 행렬 공식의 값에 따라 고정점의 안정성을 분류한 도표이다.

게임에서 각양각색의 2변수의 비선형 역학 시스템의 수치를 계산할 때, 각 방정식에 대해 해석적으로 미분함수 $f (x, y), g (x, y) $를 편향시켜 설명하는 것은 좀 번거롭다.이에 따라 이번에는 프로그램이 편미분 계산을 하기 위해 쌍수 자동 미분을 채택했다.python에서 이중수를 사용하기 위해 만든 모듈 (dual.py) 의 원본 파일은 여기. 에 두어야 합니다.다음은 임의의 $x, y$의 대각과 행렬 공식을 계산하기 위한 함수입니다.

from dual import *

def calc_trace_determinant2(x,y):
    det = fx(x,y)*gy(x,y)-fy(x,y)*gx(x,y)
    tr = fx(x,y)+gy(x,y)
    return tr,det

def fx(x,y):
    return dual(f(dual(x,1),y)-f(x,y)).im

def fy(x,y):
    return dual(f(x,dual(y,1))-f(x,y)).im

def gx(x,y):
    return dual(g(dual(x,1),y)-g(x,y)).im

def gy(x,y):
    return dual(g(x,dual(y,1))-g(x,y)).im

# Brusselator 方程式の場合, f, g は次のように定義する.
# 引数x, yにはスカラーを指定してください.
def f(x,y):
    return a-(b+1)*x+x**2*y

def g(x,y):
    return b*x-x**2*y

그 밖에 선형 안정성 분석에서 고정점의 대각과 행렬 공식의 값을 알기 위해 뉴턴법으로 아래 함수를 통해 고정점의 좌표를 얻었다.

import numpy as np

def newton_method2(f,g,x_range,y_range,epsilon=10**(-13),limit=10**4):
    x_min,x_max = x_range
    y_min,y_max = y_range
    x = np.random.rand()*(x_max-x_min)+x_min
    y = np.random.rand()*(y_max-y_min)+y_min
    err_x,err_y = 1,1
    count = 0
    while err_x>epsilon or err_y>epsilon:
        _,det = calc_trace_determinant2(x,y)
        dx = (fy(x,y)*g(x,y)-f(x,y)*gy(x,y))/det
        dy = (gx(x,y)*f(x,y)-g(x,y)*fx(x,y))/det
        x = x+dx
        y = y+dy
        err_x = np.abs(dx)
        err_y = np.abs(dy)
        count += 1
        if count>limit:
            raise RuntimeError("count exceeded limit @ newton_method2.")
    return x,y

다음 그림은 상기 함수를 사용하여 Brusselator 방정식의 고정점 아크비의 대각과 행렬 공식을 계산하고 인자 $a, b$$$$$$$$$$$$tr-det 평면을 바꾸는 데 그려집니다.

여기에 위의 그림에서 달러는 아래 공식에 따라 계산된 것이다.

\begin{align}
&\hat{a} = \frac{a-a_{min}}{a_{max}-a_{min}},\qquad
\left\{\begin{aligned}
&a\in A=\{a_m\}_{1\le m\le M}\\
&a_{min}=\mathrm{min}\;A\\
&a_{max}=\mathrm{max}\;A
\end{aligned}\right.\\
&\hat{b} = \frac{b-b_{min}}{b_{max}-b_{min}},\qquad
\left\{\begin{aligned}
&b\in B=\{b_n\}_{1\le n\le N}\\
&b_{min}=\mathrm{min}\;B\\
&b_{max}=\mathrm{max}\;B
\end{aligned}\right.\\
&color = \frac{M\hat{a}+\hat{b}}{M+1},\qquad 0\le color\le 1
\end{align}

수치 계산에서 $a{min}=b_{min}=0.1$, $a_{max}=b_{max]=2달러, $M=N=10달러입니다.
tr-det 평면만 그리면 자동 미분의 은혜를 느끼지 못하기 때문에 마지막으로 x-y 평면에 고정점을 그려 매개 변수가 변할 때 고정점의 안정성이 어떤gif 애니메이션으로 변할지 색깔로 표시했다.gif 애니메이션에서 고정점의 안정성은 검은색 안정성, 회색 중심, 흰색 불안정이다.

실제 온라인 안정성 분석에서dual.py를 사용할 때 sample2.py를 참고하십시오.