WKB近似

WKB方法用于研究一种特定类型的微分方程的全局性质
- 很有用
- 这种特定的微分方程形如：

$y''+[\lambda^2q_1(x)+q_2(x)]y=0$

经过一些不是特别复杂的推导，我们可以得到他的WKB近似解。
- 该近似解的选择取决于函数 $q_1(x)$ 和参数 $\lambda$ 的性质
- 同时，我们默认函数的定义域为 $[0,+\infty)$
当 $q_1(x)$ 恒大于零, $\lambda\rightarrow +\infty$ 时：

$y=\frac{1}{q_1^{\frac{1}{4}}(x)}\begin{pmatrix} a\cdot cos(\lambda \int_0^x\sqrt{q_1(\tau)}d\tau)+b\cdot sin(\lambda\int_0^x\sqrt{q_1(\tau)}d\tau) \end{pmatrix}$

当 $q_1(x)$ 恒小于零, $\lambda\rightarrow +\infty$ 时：

$y=\frac{1}{[-q_1(x)]^{\frac{1}{4}}} \begin{pmatrix} a\cdot ch(\lambda \int_0^x\sqrt{-q_1(\tau)}d\tau)+b\cdot sh(\lambda\int_0^x\sqrt{-q_1(\tau)}d\tau) \end{pmatrix}$

注：很快你会发现两个问题
- a 和 b 是怎么确定的
  - 答：初始条件
- \lambda 趋于正无穷并不意味着只有 \lambda 很大才能使得该WKB 近似真实
当存在一个或者多个零点时，函数的WKB近似会变得非常的复杂
- 我们一会再看

WKB 近似案例 1

该案例选取自参考文献【1】P137 例5.1.1

$\begin{matrix} y''=\lambda^2Q(x)y\\ y(0)=0,y'(0)=1\\ Q(x)=(1+x^2)^2 \end{matrix}$

我们相信它的WKB近似为 $y\sim \frac{1}{\lambda (1+x^2)^{\frac{1}{2}}}sh(\lambda\cdot(x+\frac{x^3}{3}))$
我们使用四阶Runge-Kutta方法求解

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fsolveclass odessolver():def __init__(self, f, Y_start=np.array([0, 1]), dY_start=np.array([0, 0]), \X_start=0, X_end=1, h=0.01):self.f = fself.h = hself.X = np.arange(X_start, X_end, self.h)self.n = Y_start.sizeself.Y = np.zeros((self.n, self.X.size))#第一个参数表示元 第二个参数表示变量self.Y[:, 0] = Y_startself.Y[:, 1] = Y_start + self.h * dY_startself.tol = 1e-6def __str__(self):return f"y'(x) = f(x) = ({self.f}) variables"def RK4(self):for i in range(1, self.X.size):k1 = self.f(self.X[i-1]           , self.Y[:, i-1])k2 = self.f(self.X[i-1] +self.h/2 , self.Y[:, i-1]+1/2*self.h*k1)k3 = self.f(self.X[i-1] +self.h/2 , self.Y[:, i-1]+1/2*self.h*k2)k4 = self.f(self.X[i-1] +self.h   , self.Y[:, i-1]+    self.h*k3)self.Y[:, i] = self.Y[:, i-1] +self.h/6 * (k1 + 2*k2 + 2*k3 + k4)return self.Ydef IRK4(self):for i in range(1, self.X.size):def f1(k1, k2):f1_x = self.X[i-1] + self.h*(3-3**0.5)/6f1_y = self.Y[:, i-1]+k1/4*self.h+(3-2*3**0.5)/12*k2*self.hf1_res = self.f(f1_x, f1_y)return np.array([f1_res[i] for i in range(self.n)])def f2(k1, k2):f2_x = self.X[i-1] + self.h*(3+3**0.5)/6f2_y = self.Y[:, i-1]+k2/4*self.h+(3+2*3**0.5)/12*k1*self.hf2_res = self.f(f2_x, f2_y)return np.array([f2_res[i] for i in range(self.n)])              def func(k):k1 = np.array([k[i] for i in range(self.n)])k2 = np.array([k[i+self.n] for i in range(self.n)])doc = []for i in range(self.n):doc.append((k1 - f1(k1, k2))[i])for i in range(self.n):doc.append((k2 - f2(k1, k2))[i])return docsol = fsolve(func, np.zeros(self.n*2))self.Y[:, i] = self.Y[:, i-1] + 1/2 * self.h * (sol[:self.n] + sol[self.n:])return self.YA = 0 
B = 1
Lambda = 1
Q = lambda x:(1+x**2)**2
Y0 = np.array([A, B])def test_fun(x, Y):   return np.array([Y[1], Lambda**2 * Q(x) * Y[0]])c = odessolver(test_fun, Y_start=Y0)
x = np.arange(0, 1, 0.01)y3 = c.RK4()
x = np.arange(0, 1, 0.01)
plt.plot(x, y3[0, :], label="RK4")##y4 = c.IRK4()
##x = np.arange(0, 1, 0.01)
##plt.plot(x, y4[0, :], label="IRK4")WKB = lambda x:1/(Lambda*(1+x**2)**0.5)*(np.exp(x+x**3/3)-np.exp(-(x+x**3/3)))/2
plt.plot(x, WKB(x), label="WKB")plt.legend()
plt.pause(0.01)