title Sturm-Liouville theory jupyter python3 We cn use the `scipy.integrate.solve_bvp` method to solve many (most?) boundary value problems, including eigenvalue problems (in which an additional parameter must be determined). ```{python} import numpy as np from scipy.integrate import solve_bvp import scipy as sp import matplotlib.pyplot as plt ``` ## Sheet 1, question 6 The system is: $$ y'' + \lambda y = 0, \quad y(0) = 0, \quad y(1) + y'(1) = 0 $$ We can write this system in the form $$ (-D^2) y = \lambda y $$ where the differential operator $D = \mathrm{d}/\mathrm{d}x$ will be treated as though it were a matrix. ```{python} def fun(x, y, p): lamb = p[0] return np.vstack((y[1], -lamb * y[0])) def bc(ya, yb, p): lamb = p[0] return np.array([ya[0], yb[0] + yb[1], ya[1] - 1]) mode = 2 x_guess = np.linspace(0, 1, 101) y_guess = np.zeros((2, x_guess.size)) y_guess[0] = np.sin(mode*np.pi*x_guess) / (mode*np.pi) lamb_guess = (mode*np.pi) ** 2 sol = solve_bvp(fun, bc, x_guess, y_guess, p=[lamb_guess]) print(sol.p[0]) x_plot = np.linspace(0, 1, 100) y_plot = sol.sol(x_plot)[0] plt.plot(x_plot, y_plot) plt.title(f"lambda = {sol.p[0]}") plt.grid() plt.xlabel("x") plt.ylabel("y") plt.show() ``` ## Sheet 1, question 7 $$ y'' + 4y' + (4+\lambda)y = 0, \quad y(0) = y(1) = 0 $$ ```{python} def fun(x, y, p): lamb = p[0] v = (y[1], -y[1] - (4+lamb)*y[0]) return np.vstack(v) def bc(ya, yb, p): lamb = p[0] # y'(0) = 1 normalisation return np.array([ya[0], yb[0], ya[1] - 1]) # The mode that gets found depends heavily on the initial guess, # which you can use to set the approximate location of the first # crest. mode = 3 lamb_guess = 0 x_guess = np.linspace(0, 1, 101) y_guess = np.zeros((2, x_guess.size)) y_guess[0, :] = np.sin(mode * x_guess * np.pi) / (mode*np.pi) sol = solve_bvp(fun, bc, x_guess, y_guess, p=[lamb_guess]) x_plot = np.linspace(0, 1, 100) y_plot = sol.sol(x_plot)[0] plt.plot(x_plot, y_plot, 'k-', x_guess, y_guess[0], 'k:') plt.title(f"lambda = {sol.p[0]}") plt.legend(["solution", "initial guess"]) plt.grid() plt.xlabel("x") plt.ylabel("y") plt.show() ``` ## Sheet 1, question 8: Bessel's equation (RSP at $x=0$) Bessel's equation arises when solving Laplace's equation in cylindrical coordinates using the method of separation. It determines the radial behaviour of the solution. The eigenvalues $\lambda$ come from the contributions of the $\theta$ and $z$ terms in the Laplacian. $$ (x y')' + \lambda x y = 0, \quad y(0) = 1, y(1) = 0 $$ Specifically, this is Bessel's equation _of order zero_, which arises for a cylindrically symmetric solution. Recast the LHS as a singular DE: $$ y'' + \frac{1}{x} y' + \lambda y = 0 $$ For normalisation we'll take $y(0) = 1$. This is very similar to the harmonic equation $y'' + \lambda y = 0$; the $y'/x$ term from the curvilinear coordinate system gets less and less important as $x\rightarrow\infty$. ```{python} def fun(x, y, p): lamb = p[0] v = (y[1], y[2], -y[1]/(x+0.01) - lamb*y[0]) return np.vstack(v) def bc(ya, yb, p): lamb = p[0] # y(0) = 1 normalisation return np.array([ya[0] - 1, ya[1], ya[2] - 1, yb[0]]) # The mode that gets found depends heavily on the initial guess, # which you can use to set the approximate location of the first # crest. mode = 2 lamb_guess = 15 x_guess = np.linspace(0, 1, 101) y_guess = np.zeros((3, x_guess.size)) y_guess[0, :] = np.cos((mode+1/2) * x_guess * np.pi) sol = solve_bvp(fun, bc, x_guess, y_guess, p=[lamb_guess]) x_plot = np.linspace(0, 1, 100) y_plot = sol.sol(x_plot)[0] plt.plot(x_plot, y_plot, 'k-', x_guess, y_guess[0], 'b:') plt.title(f"lambda = {sol.p[0]}") plt.legend(["solution", "initial guess"]) plt.grid() plt.xlabel("x") plt.ylabel("y") plt.show() ``` ## Question 9: Higher-order self-adjoint form $$ - y'''' + \lambda y = 0 $$ with homogeneous BCs in $y, y'$. ```{python} def fun(x, y, p): lamb = p[0] return np.vstack((y[1], y[2], y[3], lamb * y[0])) def bc(ya, yb, p): lamb = p[0] return np.array([ya[0], ya[1], yb[0], yb[1], ya[2] - 1]) mode = 4 x_guess = np.linspace(0, 1, 101) y_guess = np.zeros((4, x_guess.size)) y_guess[0, :] = x_guess * (1-x_guess) * np.sin(mode * x_guess * np.pi) / (2*mode*np.pi) sol = solve_bvp(fun, bc, x_guess, y_guess, p=[12]) x_plot = np.linspace(0, 1, 100) y_plot = sol.sol(x_plot)[0] # Normalise the guess (just for demonstration) plt.plot(x_plot, y_plot, 'k-', x_guess, y_guess[0] * np.max(np.abs(y_plot)) / np.max(np.abs(y_guess)), 'k:') plt.title(f"lambda = {sol.p[0]}") plt.grid() plt.xlabel("x") plt.ylabel("y") plt.show() ```