\[\frac{d^{2}}{d r^{2}} \phi{\left(r \right)} = - 0.5 r e^{- r}\]
Chapt 3: BVP & EVP of ODE
Feng Wan 
B832@仲英楼 or 3133@泓理楼 // School of Physcis @ Xian Jiaoton University
Code
def ode(x, y):
res = y.copy()
res[0] = y[1]
res[1] = -4*np.pi**2 * y[0]
return res
def rk4(func, x0, y0, h, xp):
steps = int((xp - x0) / h)
result = np.zeros((steps + 1, len(y0)))
xis = np.zeros((steps + 1, 1))
xi = x0
yi = y0
xis[0] = x0
result[0, :] = yi
for i in range(steps):
k1 = func(xi, yi)
k2 = func(xi + 0.5 * h, yi + 0.5 * h * k1)
k3 = func(xi + 0.5 * h, yi + 0.5 * h * k2)
k4 = func(xi + h, yi + h * k3)
xi = xi + h
yi = yi + h / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
xis[i + 1] = xi
result[i + 1, :] = yi
return xis, result
h = 0.01
x0 = 0
xp = 1.0
y0 = np.array([1.0, 0.0])
xi, yi = rk4(ode, x0, y0, h, xp)
plt.plot(xi, yi[:, 0], label = r'$y$')
plt.plot(xi, yi[:, 1], label = r'$dy/dt$')
plt.xlabel(r'$t$')
plt.legend()
plt.show()
Code
def numerov(x0, y0, y0p, h, xp):
steps = int((xp - x0) / h)
result = np.zeros((steps + 1, 1))
xis = np.zeros((steps + 1, 1))
y1 = y0 + h * y0p + (- 4.0 * np.pi**2.0 * y0) * h**2.0 / 2.0
result[1] = y1
result[0] = y0
fac = 2.0 * (3.0 - 5.0 * np.pi**2.0 * h**2.0) / (3.0 + np.pi**2.0 * h**2.0)
xis[0] = x0
xis[1] = x0 + h
for i in range(steps - 1):
j = i + 1
xis[j + 1] = (j + 1) * h + x0
result[j + 1] = fac * result[j] - result[j - 1]
return xis, result
def numerov2(x0, y0, y1, h, xp):
steps = int((xp - x0) / h)
result = np.zeros((steps + 1, 1))
xis = np.zeros((steps + 1, 1))
result[1] = y1
result[0] = y0
fac = 2.0 * (3.0 - 5.0 * np.pi**2.0 * h**2.0) / (3.0 + np.pi**2.0 * h**2.0)
xis[0] = x0
xis[1] = x0 + h
for i in range(steps - 1):
j = i + 1
xis[j + 1] = (j + 1) * h + x0
result[j + 1] = fac * result[j] - result[j - 1]
return xis, result
x0 = 0
y0 = 1
y0p = 0
h = 0.01
xp = 1
xis, result = numerov(x0, y0, y0p, h, xp)
xis2, result2 = numerov2(x0, y0, yi[1, 0], h, xp)
ex_y = np.cos(2.0 * np.pi * xi[:, 0])
plt.plot(xi, ex_y, label = 'explicit')
plt.plot(xis, result, 'r--', label = 'numerov')
plt.plot(xi, yi[:, 0], 'b-.', label = 'rk4')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()
Code
def ode1(x, y):
res = y.copy()
res[0] = y[1]
res[1] = -x * y[1] + 4 * y[0] + 12 * x**2 - 3 * x
return res
def ode2(x, y):
res = y.copy()
res[0] = y[1]
res[1] = - x * y[1] + 4 * y[0]
return res
y00 = np.array([0, 0])
y01 = np.array([0, 1])
h = 0.02
x0 = 0
xp = 1
xi1, yi1 = rk4(ode1, x0, y00, h, xp)
xi2, yi2 = rk4(ode2, x0, y01, h, xp)
y1 = yi1[:, 0]
y2 = yi2[:, 0]
yfin = y1 + (2.0 - y1[-1]) / y2[-1] * y2
y_ex = xi1**4.0 + xi1
plt.plot(xi1, yfin, 'r-', label = 'superstition')
plt.plot(xi1, y_ex, 'b-.', label = 'explicit')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()
Code
def shooting(ode, x0, y0, h, xp, yp, eps):
s1 = 0.1
s2 = 0.2
y00 = np.array([y0, 0.0])
while true:
y00[1] = s1
x1, y1 = rk4(ode, x0, y00, h, xp)
y00[1] = s2
x2, y2 = rk4(ode, x0, y00, h, xp)
y11 = y1[:, 0]
y22 = y2[:, 0]
if (np.abs(y11[-1] - yp) < eps):
return x1, y11
if (np.abs(y22[-1] - yp) < eps):
return x2, y22
s3 = s2 - (s2 - s1) / (y22[-1] - y11[-1]) * (y22[-1] - yp)
print(s3)
y00[1] = s3
x3, y3 = rk4(ode, x0, y00, h, xp)
y33 = y3[:, 0]
if (np.abs(y33[-1] - yp) < eps):
return x3, y33
s1 = s2
s2 = s3
def ode(x, y):
res = y.copy()
res[0] = y[1]
res[1] = -x * y[1] + 4 * y[0] + 12 * x**2 - 3 * x
return res
x0 = 0
xp = 1
y0 = 0
yp = 2
h = 0.02
eps = 1e-6
xx, yy = shooting(ode, x0, y0, h, xp, yp, eps)
plt.plot(xx, yy)
plt.xlabel('x')
plt.ylabel('y')
plt.show()1.0000000334999262
Code
def shooting(ode, x0, y0, h, xp, yp, eps):
s1 = 0.1
s2 = 0.2
y00 = np.array([0.0, y0])
while true:
y00[0] = s1
x1, y1 = rk4(ode, x0, y00, h, xp)
y00[0] = s2
x2, y2 = rk4(ode, x0, y00, h, xp)
y11 = y1[:, 0]
y22 = y2[:, 0]
if (np.abs(y1[-1, 1] - yp) < eps):
return x1, y11
if (np.abs(y2[-1, 1] - yp) < eps):
return x2, y22
s3 = s2 - (s2 - s1) / (y2[-1, 1] - y1[-1, 1]) * (y2[-1, 1] - yp)
y00[0] = s3
x3, y3 = rk4(ode, x0, y00, h, xp)
y33 = y3[:, 0]
if (np.abs(y3[-1, 1] - yp) < eps):
return x3, y33
s1 = s2
s2 = s3
def ode(x, y):
res = y.copy()
res[0] = y[1]
res[1] = -x * y[1] + 4 * y[0] + 12 * x**2 - 3 * x
return res
x0 = 0
xp = 1.0
y0 = 1.0
yp = 5.0
h = 0.02
eps = 1e-6
xx, yy = shooting(ode, x0, y0, h, xp, yp, eps)
plt.plot(xx, yy)
plt.xlabel('x')
plt.ylabel('y')
plt.show()
分离变量后，空间部分满足的方程和边界条件可以写成 \[ \left\lbrace \begin{aligned}[l] &\partial_x^2 \phi = -k^2\phi; \\ &\phi(0) = \phi(1) = 0 \end{aligned} \right. \] \(k\) wavenumber, \(\phi\) transverse motion of the string.
Explicit solution: \[
\begin{aligned}[l]
&k_n = n\pi; \\
&\phi_n \sim \sin n\pi x;~~ (n = 1, 2, ...)
\end{aligned}
\] 相比边值问题，本征值问题多了一个待定参数 ——
Code
eps = 1e-4
x0 = 0.0
xp = 1.0
h = 0.01
y0 = 0.0
yp = 0.0
def ode(x, y, k):
res = y.copy()
res[0] = y[1]
res[1] = -k*k*y[0]
return res
def shooting_evp(ode, x0, y0, h, xp, yp, eps):
k1 = 0.5
k2 = 5.0
delta = 30.0
y00 = np.array([y0, delta])
node = lambda x, y: ode(x, y, k1)
x1, y1 = rk4(node, x0, y00, h, xp)
f1 = y1[-1, 0] - yp
node = lambda x, y: ode(x, y, k2)
x2, y2 = rk4(node, x0, y00, h, xp)
f2 = y2[-1, 0] - yp
if (f1 * f2 > 0):
print("no roots in the interval")
exit(1)
if (f1 * f2 < 0):
if (np.abs(k2 - k1) < eps):
return 0.5 * (k1 + k2), xx, yy
while True:
kk = 0.5 * (k1 + k2)
node = lambda x, y: ode(x, y, kk)
xx, yy = rk4(node, x0, y00, h, xp)
plt.plot(xx, yy[:, 0], '-.')
if (np.abs(k2 - k1) < eps):
return kk, xx, yy
fx = yy[-1, 0] - yp
if (f1 * fx < 0):
k2 = kk
else:
k1 = kk
node = lambda x, y: ode(x, y, k1)
x1, y1 = rk4(node, x0, y00, h, xp)
f1 = y1[-1, 0] - yp
kval, xval, yval = shooting_evp(ode, x0, y0, h, xp, yp, eps)
plt.grid('on')
plt.plot(xval, yval[:, 0], 'k-')
plt.show()
Code
eps = 1e-4
x0 = 0.0
xp = 1.0
h = 0.1
y0 = 0.0
yp = 0.0
def ode(x, y, k):
res = y.copy()
res[0] = y[1]
res[1] = -k*k*y[0]
return res
def shooting_evp(ode, x0, y0, h, xp, yp, eps, n):
j = 1
k1 = 0.0
k2 = k1 + h
delta = 30.0
y00 = np.array([y0, delta])
node = lambda x, y: ode(x, y, k1)
x1, y1 = rk4(node, x0, y00, h, xp)
f1 = y1[-1, 0] - yp
node = lambda x, y: ode(x, y, k2)
x2, y2 = rk4(node, x0, y00, h, xp)
f2 = y2[-1, 0] - yp
if (f1 * f2 > 0):
k1 = k2
k2 = k1 + h
if (f1 * f2 < 0):
if (np.abs(k2 - k1) < eps):
print(j, kk)
plt.plot(xx, yy[:, 0], label = 'k = ' + str(j))
j = j + 1
while j < n:
kk = 0.5 * (k1 + k2)
node = lambda x, y: ode(x, y, kk)
xx, yy = rk4(node, x0, y00, h, xp)
plt.plot(xx, yy[:, 0], '-.')
if (np.abs(k2 - k1) < eps):
print(j, kk)
plt.plot(xx, yy[:, 0], label = 'k = ' + str(j))
j = j + 1
k1 = k2
k2 = k1 + h
fx = yy[-1, 0] - yp
if (f1 * fx < 0):
k2 = kk
else:
k1 = kk
node = lambda x, y: ode(x, y, k1)
x1, y1 = rk4(node, x0, y00, h, xp)
f1 = y1[-1, 0] - yp
shooting_evp(ode, x0, y0, h, xp, yp, eps, 3)
plt.grid('on')
plt.legend()
plt.show()1 0.199951171875
2 0.29995114803314216
