Inner Product
If \(\,\mathbf{u}=u_1\mathbf{i} +u_2\mathbf{j} +u_3\mathbf{k}\,\) and \(\,\mathbf{v} =v_1\mathbf{i} +v_2\mathbf{j} +v_3\mathbf{k}\,\) are two vectors in \(\mathbb{R}^3\), \(\,\)then inner product is defined:
\[(\mathbf{u},\mathbf{v})=u_1v_1 +u_2v_2 +u_3v_3=u_iv_i\]
The inner product \((\mathbf{u},\mathbf{v})\) possesses the following properties:
\((\mathbf{u},\mathbf{v})=(\mathbf{v},\mathbf{u})\)
\((k\mathbf{u},\mathbf{v}) =k(\mathbf{u},\mathbf{v})\), \(\;k:\;\) a scalar
\((\mathbf{u},\mathbf{u})=0\,\) if \(\,\mathbf{u}=\mathbf{0}\,\) and \(\,(\mathbf{u},\mathbf{u})>0\,\) if \(\,\mathbf{u}\neq\mathbf{0}\)
\((\mathbf{u}+\mathbf{v},\mathbf{w})=(\mathbf{u},\mathbf{w}) +(\mathbf{v},\mathbf{w})\)
Inner Product of Functions
The inner product of two functions \(f_1\) and \(f_2\) on an interval \([a,b]\) is the number
\[ (f_1,f_2)=\int_a^b f_1(x)\, f_2(x) \,dx\]
Orthogonal Functions
Two functions \(\,f_1\) and \(\,f_2\) are said to be orthogonal on \([a,b]\,\) if
\[(f_1,f_2)=\int_a^b f_1(x)\, f_2(x) \,dx=0\]
Orthogonal Sets
A set of real-valued functions \(\{\phi_0(x),\phi_1(x),\phi_2(x),\cdots\}\,\) is said to be orthogonal on \([a,b]\,\) if
\[(\phi_m,\phi_n)=\int_a^b \phi_m(x) \, \phi_n(x) \,dx=0, \;\;m\neq n\]
With the property that \(\|\phi_n(x)\|=1\), \(\,\)then \(\{\phi_n(x)\}\) is said to be an orthonormal set on the interval
\(~\)
Example \(\displaystyle \left\{ \frac{1}{\sqrt{2\pi}}, \frac{\cos x}{\sqrt{\pi}}, \frac{\cos 2x}{\sqrt{\pi}},\cdots \right\}\), \(\;\; -\pi \leq x\leq \pi\)
\(~\)
Orthogonal Series Expansion
Suppose \(\left\{ \phi_n(x)\right\}\) is an infinite orthogonal set of functions on \([a,b]\). \(\,\)If \(\,y=f(x)\,\) is a function defined on \([a,b]\),
\[ \begin{aligned} \color{blue}{f(x)} &\; \color{blue}{= c_0\phi_0(x) +c_1\phi_1(x) +\cdots +c_n\phi_n(x)+\cdots}\\ &\; {\scriptsize\Downarrow}\\ {\scriptsize \int_a^b f(x) \phi_n(x)\,dx} &\; {\scriptsize =c_0\int_a^b \phi_0(x) \phi_n(x)\,dx +c_1\int_a^b \phi_1(x)\phi_n(x)\,dx +\cdots +c_n\int_a^b \phi_n^2(x)\,dx+\cdots}\\ &\; {\scriptsize = c_n\int_a^b \phi_n^2(x)\,dx} \\ &\; {\scriptsize\Downarrow} \\ {\color{red}{c_n}} & {\color{red}{= \frac{\displaystyle \int_a^b f(x)\,\phi_n(x)\,dx}{\displaystyle \int_a^b \phi_n^2(x)\,dx},\;\;n=0,1,2,\cdots}} \end{aligned}\]
Orthogonal Set/Weight Function
A set of real-valued functions \(\{\,\phi_0(x),\phi_1(x),\phi_2(x),\cdots\}\,\) is said to be orthogonal with respect to a weight function \(w(x)\) on \([a,b]\,\) if
\[(\phi_m,\phi_n)=\int_a^b w(x)\phi_m(x) \phi_n(x) \,dx=0, \;\;m\neq n\]
With this inner product notation,
\[ f(x)=\sum_{n=0}^\infty c_n \phi_n(x)\]
where \(~{\scriptsize c_n=\frac{\displaystyle\int_a^b f(x)w(x)\phi_n(x)\,dx}{\displaystyle\int_a^b w(x) \phi_n^2(x)\,dx}}\)
\(~\)
Trigonometric Series
\[ \left\{ 1, \,\cos\frac{\pi}{p}x, \,\cos\frac{2\pi}{p}x, \,\cos\frac{3\pi}{p}x,\cdots, \,\sin\frac{\pi}{p}x, \,\sin\frac{2\pi}{p}x, \,\sin\frac{3\pi}{p}x, \,\cdots \right\} \]
Fourier Series
The Fourier series of a function \(f\) defined on the interval \((-p,p)\,\) is given by
\[ f(x)=\frac{a_0}{2} +\sum_{n=1}^\infty \left( a_n\cos\frac{n\pi}{p}x +b_n\sin\frac{n\pi}{p}x \right)\]
where
\[ \begin{aligned} a_n &= \frac{1}{p}\int_{-p}^p f(x)\cos\frac{n\pi}{p}x \,dx, \quad a_0 = \frac{1}{p}\int_{-p}^p f(x)\,dx \\ b_n &= \frac{1}{p}\int_{-p}^p f(x)\sin\frac{n\pi}{p}x \,dx \end{aligned}\]
At a point of discontinuity, \(\,\)the Fourier series converges to the average \(\displaystyle {\scriptsize \;\frac{f(x^+)+f(x^-)}{2}}\)
\(~\)
Example
\[\begin{aligned} f(x) &= \begin{cases} \;0, & -\pi \leq x \leq 0 \\ \;1, & \;\;\, 0 < x \leq \pi \end{cases} \\ &\Downarrow \\ f(x)&= \frac{1}{2} +\sum_{m=1}^\infty \frac{2}{(2m -1)\pi} \sin{(2m -1)x} \end{aligned}\]
\(~\)
from sympy import fourier_series, Piecewise, pi, plot
from sympy.abc import x
f = Piecewise((0, x < 0), (1, True))
s = fourier_series(f, (x, -pi, pi))
s.truncate(5)\(\displaystyle \frac{2 \sin{\left(x \right)}}{\pi} + \frac{2 \sin{\left(3 x \right)}}{3 \pi} + \frac{2 \sin{\left(5 x \right)}}{5 \pi} + \frac{2 \sin{\left(7 x \right)}}{7 \pi} + \frac{1}{2}\)
s3 = s.truncate(n = 3)
s6 = s.truncate(n = 6)
s9 = s.truncate(n = 9)
p = plot(f, s3, s6, s9, (x, -1.5 *pi, 1.5 *pi), show=False, legend=True)
p[0].label = 'f'
p[1].label = 'n = 3'
p[2].label = 'n = 6'
p[3].label = 'n = 9'
p.show()
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 10
x_t1 = [-2*np.pi, -np.pi, -np.pi, 0, 0, np.pi, np.pi, 2*np.pi]
y_t1 = [1, 1, 0, 0, 1, 1, 0, 0]
x_t2 = [-np.pi, 0, 0, np.pi]
y_t2 = [0, 0, 1, 1]
x = np.linspace(-2 *np.pi, 2 *np.pi, 1000)
def example_plot(func, ns):
fig, axes = plt.subplots(2, 2, figsize=(7.5, 7))
i = 0
for n in ns:
ax = axes[i//2, i%2]
ax.plot(x_t1, y_t1, 'b:')
ax.plot(x_t2, y_t2, 'g-', lw=3)
ax.plot(x, func(x, n), c='darkorange')
ax.set_xticks([-2.0 *np.pi, -np.pi, 0, np.pi, 2.0 *np.pi])
ax.set_xticklabels(['$-2\pi$', '$-\pi$', '0', '$\pi$', '$2\pi$'])
ax.set_xlim(-2.0 *np.pi, 2.0 *np.pi)
ax.set_ylim(-0.2, 1.2)
ax.set_title(f"$n=${n}")
i += 1def rect_wave(x, n):
y = 0.5
for m in range(1, n +1):
y = y +2.0 /((2.0*m -1) *np.pi) *np.sin((2.0*m -1)*x)
return y
example_plot(rect_wave, ns=[1, 3, 6, 9])
\(~\)
\(~\)
If \(\,f\,\) is even, \(~\) then
\(\displaystyle\phantom{xx} \int_{-a}^a f(x)\,dx = 2\int_0^a f(x)\,dx\)
If \(\,f\,\) is odd, \(\,\)then
\(\displaystyle\phantom{xx} \int_{-a}^a f(x)\,dx = 0\)
The Fourier series of an even function on the interval \((-p,p)\,\) is the cosine series
\[ f(x)=\frac{a_0}{2} +\sum_{n=1}^\infty a_n \cos\frac{n\pi}{p}x \]
\(\text{where }\,\displaystyle a_0=\frac{2}{p} \int_0^p f(x)\,dx,\; a_n=\frac{2}{p} \int_0^p f(x) \cos\frac{n\pi}{p}x \,dx\)
The Fourier series of an odd function on the interval \((-p,p)\,\) is the sine series
\[ f(x)=\sum_{n=1}^\infty b_n \sin\frac{n\pi}{p}x\;\;\;\text{where}\;\displaystyle b_n=\frac{2}{p} \int_0^p f(x) \sin\frac{n\pi}{p}x \,dx \]
\(~\)
Example
\[ \begin{aligned} f(x) &= x, \;\;0 < x < 2 \\ &\Downarrow \\ f(x) &= \frac{4}{\pi}\sum_{n=1}^\infty \frac{(-1)^{n +1}}{n} \sin{\frac{n\pi}{2}x} \\ &\,\text{or} \\ f(x) &= 1 -\frac{8}{\pi^2}\sum_{m=1}^\infty \frac{1}{(2m -1)^2} \cos{\frac{(2m -1)\pi}{2}x} \end{aligned}\]
\(~\)
x_t1 = [-4,-2,-2, 2, 2, 4]
y_t1 = [ 0, 2,-2, 2,-2, 0]
x_t2 = [0, 2]
y_t2 = [0, 2]
x = np.linspace(-4, 4, 1000)
def example_plot2(func, ns):
fig, axes = plt.subplots(2, 2, figsize=(7.5, 7))
i = 0
for n in ns:
ax = axes[i//2, i%2]
ax.plot(x_t1, y_t1, 'b:')
ax.plot(x_t2, y_t2, 'g-', lw=3)
ax.plot(x, func(x, n), c='darkorange')
ax.set_xlim(-4, 4)
ax.set_title(f"$n=${n}")
i += 1def lamp_sin(x, n):
y = 0
for m in range(1, n +1):
y = y +4.0 /np.pi *(-1)**(m +1)/m *np.sin(m *np.pi /2 *x)
return y
example_plot2(lamp_sin, ns=[1, 3, 6, 9])
x_t1 = [-4,-2, 0, 2, 4]
y_t1 = [ 0, 2, 0, 2, 0]
def lamp_cos(x, m):
y = 1
for m in range(1, m +1):
y = y -8.0 /(np.pi *(2*m -1))**2.0 *np.cos((2*m -1) *np.pi /2 *x)
return y
example_plot2(lamp_cos, ns=[1, 3, 6, 9])
\(~\)
Half-Range Expansion
Throughout the preceding discussion, \(\,\)it was understood that a function \(\,f\,\) was defined on an interval with the origin as midpoint
However, in many instances, \(\,\)we are interested in representing a function that is defined only for \(\,0 < x < L~\) by a trigonometric series
This can be done in many different ways by supplying an arbitrary definition of the function on the interval \(-L < x < 0\). \(\,\)For brevity we condider three most important cases:
Reflect the graph of the function about the \(y\)-axis onto \(-L < x < 0\)
Reflect the graph of the function through the origin onto \(-L < x < 0\)
Define \(f\) on \(-L < x < 0\;\) by \(\,f(x)=f(x+L)\)
Complex Fourier Series
The complex Fourier series of functions \(\,f\) defined on an interval \((-p,p)\,\) is given by
\[ f(x)=\sum_{n=-\infty}^\infty c_n e^{\frac{in\pi}{p}x} \]
where
\[ c_n=\frac{1}{2p}\int_{-p}^p f(x) e^{-\frac{in\pi}{p}x}\,dx \]
Fundamental Frequency
The Fourier series define a periodic function and the fundamental period of that function is \(T=2p\). \(\,\)Since \(p=T/2\),
\[ \frac{a_0}{2} +\sum_{n=1}^\infty (a_n \cos n\omega x +b_n \sin n\omega x)\;\;\]
\[\text{and}\]
\[ \sum_{n=-\infty}^\infty c_n e^{in\omega x}\]
where \(\,\omega=2\pi/T\,\) is called the fundamental angular frequency
Frequency Spectrum
If \(\,f\) is periodic and has fundamental period \(T\), \(\,\)the plot of the points \((n\omega, |c_n|)\) \(\,\)is called the frequency spectrum of \(\;f\)
\(~\)
Example \(\,\) Expand \(\;f(x)=e^{-x}\), \(\;-\pi<x<\pi\;\) in a complex Fourier series and find the frequency spectrum
\(~\)
Self-Adjoint Form
\[\begin{aligned} a(x) \frac{d^2 y}{dx^2} +b(x) \frac{dy}{dx} &+\left[ \lambda c(x) +d(x) \right] y = 0,\;\;a(x)\neq 0 \\ &\Downarrow \\ {\scriptsize \frac{d}{dx}\left[ \exp\left(\int\frac{b(x)}{a(x)}dx\right)\frac{dy}{dx} \right ]} &{\scriptsize +\left\{ \lambda \frac{c(x)}{a(x)}\exp\left(\int\frac{b(x)}{a(x)}dx\right) +\frac{d(x)}{a(x)}\exp\left(\int\frac{b(x)}{a(x)}dx\right)\right \} y = 0} \\ &\Downarrow \\ \color{red}{\frac{d}{dx}\left[r(x)\frac{dy}{dx}\right]}&\color{red}{+\left(\lambda p(x) +q(x)\right)y=0} \end{aligned}\]
Eigenvalues and Eigenfunctions
An orthogonal set of functions can be generated by solving a two-point boundary-value problem involving a linear second-order differential equation containing a parameter \(\lambda\)
The boundary-value problem
\[ y''+\lambda y =0, \;y(0)=0, \;y(L)=0~\]
possessed nontrivial solutions only when the parameter \(\lambda\) took on the values
\[ \lambda_n= \left(\frac{n\pi}{L}\right)^2, \;\;n=1,2,3,\cdots \]
called eigenvalues. The corresponding nontrivial solutions
\[ y_n=c_n\sin \frac{n\pi x}{L} \]
are called the eigenfunctions
Regular Sturm-Liouville Problem
Let \(p\), \(q\), \(r\), and \(r'\) be real-valued functions continuous on an interval \([a,b]\), \(\,\)and let \(r(x) > 0~\) and \(\,p(x)>0~\) for every \(\,x\,\) in the interval. \(\,\)Then
\[ \frac{d}{dx}\left[r(x)\frac{dy}{dx}\right]+\left[\lambda p(x) +q(x)\right]y=0 \]
subject to
\[\begin{aligned} A_1 y(a) +B_1 y'(a)&= 0\\ A_2 y(b) +B_2 y'(b)&= 0 \end{aligned}\]
is said to be a regular Sturm-Liouville problem
Properties of the Regular Sturm-Liouville problem
There are an infinite number of real eigenvalues that can be arranged in increasing order \(\lambda_1<\lambda_2<\lambda_3<\cdots<\lambda_n<\cdots\;\;\) such that \(\lambda_n \to \infty\;\) as \(n\to\infty\)
For each eigenvalue \(\lambda_i\), \(~\)there is only one eigenfunction (except for nonzero constant multiples)
Eigenfunctions corresponding to different eigenvalues are linearly independent
The set of eigenfunctions corresponding to the set of eigenvalues is orthogonal with respect to the weight function \(p(x)\) on the interval \([a,b]\)
Proof of 4
\[\scriptstyle \begin{aligned} \frac{d}{dx}\left[r(x)\frac{dy_m}{dx}\right] &+\left(\lambda_m p(x) +q(x)\right) y_m \;=\; 0 \\ \frac{d}{dx}\left[r(x)\frac{dy_n}{dx}\right] &+\left(\lambda_n p(x) +q(x)\right) y_n \;=\; 0 \\ &\Downarrow \\ \left( \lambda_m -\lambda_n \right ) p(x) y_m y_n \; &=\; y_m \frac{d}{dx}\left[r(x)\frac{dy_n}{dx}\right] -y_n \frac{d}{dx}\left[r(x)\frac{dy_m}{dx}\right] \\ &\Downarrow \\ \left( \lambda_m -\lambda_n \right) \int_a^b p(x) y_m y_n \,dx \; & =\; \left.y_m r(x) \frac{d y_n}{dx} \right|_a^b \;-\;\int_a^b \frac{d y_m}{dx} r(x) \frac{d y_n}{dx} \,dx \;-\;\left.y_n r(x) \frac{d y_m}{dx} \right|_a^b \;+\;\int_a^b \frac{d y_n}{dx} r(x) \frac{d y_m}{dx} \,dx \\ & =\; r(b)\left[ y_m(b) y_n'(b) - y_n(b) y_m'(b) \right ] \;-\;r(a)\left[ y_m(a) y_n'(a) - y_n(a) y_m'(a) \right ] \\ & =\; 0 \\ &\Downarrow \\ \int_a^b p(x) y_m y_n \,dx \; & =\; 0, \;\;\lambda_m \neq \lambda_n \end{aligned}\]
\(~\)
Example \(\,\) Solve the boundary-value problem
\[ y'' +\lambda y=0, \;\;y(0)=0, \;\;y(1) +y'(1)=0 \]
Solution
\[ \begin{aligned} y &= c_1\cos\sqrt{\lambda}\,x +c_2\sin\sqrt{\lambda}\,x \\ \\ &\;\Downarrow \;y(0)=0,\;y(1)+y'(1)=0 \\ \\ c_1\cos 0 &+c_2\sin 0 = 0 \;\rightarrow \;c_1=0 \\ c_2\sin\sqrt{\lambda} &+c_2\sqrt{\lambda}\cos\sqrt{\lambda}=0 \;\rightarrow \;c_2 \neq 0, \;\tan\sqrt{\lambda}=-\sqrt{\lambda} \end{aligned}\]
\(~\)
x = np.linspace(0, 4*np.pi, 1000)
y = np.tan(x)
threshold = 100
y[y > threshold] = np.inf
y[y <-threshold] = np.inf
def examle_plot3(eigenvalues):
plt.figure(figsize=(6, 5))
plt.plot(x, y, linewidth=1.2, color="blue")
plt.plot(x,-x, linewidth=1.2, color="green")
plt.scatter(eigenvalues,-eigenvalues, color='red')
plt.xticks([0, np.pi, 2.0*np.pi, 3.0*np.pi, 4.0*np.pi],
['0','$\pi$','$2\pi$','$3\pi$','$4\pi$'],
fontsize=16)
plt.yticks(fontsize=16)
plt.xlim(0, 4.0*np.pi)
plt.ylim(-15, 15)
plt.grid()
plt.show()from scipy import optimize
sqrt_eigenvalues = np.zeros(4)
for i in range(4):
sqrt_eigenvalues[i] = np.round(optimize.brentq(
lambda x : np.tan(x) +x,
(i +0.5001) *np.pi,
(i +0.9999) *np.pi), 4)
print(f'{sqrt_eigenvalues =}')sqrt_eigenvalues =array([ 2.0288, 4.9132, 7.9787, 11.0855])examle_plot3(sqrt_eigenvalues)
def example_plot3_(x, n_eig, sqrt_eig, eigfunc):
plt.figure(figsize=(6, 5))
for n in range(n_eig):
plt.plot(x, eigfunc(x, sqrt_eig[n]), label='$X_%d(x)$' % (n +1))
plt.xlim(0, 1)
plt.ylim(-1.2, 1.2)
plt.grid(ls=':')
plt.legend(loc='upper left', bbox_to_anchor=(1, 1.03))
plt.xlabel('$x$')
plt.ylabel('$X_n(x)$')
plt.title(r'$X_n(x)=\sin\sqrt{\lambda_n} x$')
plt.show()
def eigfunc(x, sqrt_eig):
return np.sin(sqrt_eig *x)
x = np.linspace(0, 1, 150)
example_plot3_(x, len(sqrt_eigenvalues), sqrt_eigenvalues, eigfunc)
\(~\)
Singular Sturm-Liouville Problem
There are several important conditions under which we seek nontrivial solutions
\[\scriptsize \begin{aligned} \left( \lambda_m -\lambda_n \right ) \int_a^b p(x) y_m y_n \,dx \, & \,= r(b)\left[ y_m(b) y_n'(b) - y_n(b) y_m'(b) \right ] \\ & \text{ }-r(a)\left[ y_m(a) y_n'(a) - y_n(a) y_m'(a) \right ] = 0 \\ & \Downarrow \\ \int_a^b p(x) y_m y_n \,dx \, & \,= 0, \;\;\lambda_m \neq \lambda_n \end{aligned}\]
Parametric Bessel Equation:
\[ x^2y'' +xy' +(\color{red}{\alpha^2}x^2 -\color{blue}{n^2})y=0,\;\; n=0,1,2,\cdots\]
\[ \frac{d}{dx}\left[xy'\right] +\left(\alpha^2x -\frac{n^2}{x}\right)y=0 \]
Legendre Equation:
\[(1-x^2)y''-2xy'+\color{red}{n(n+1)}y=0,\;\; n=0,1,2,\cdots \]
\[ \frac{d}{dx}\left[(1-x^2)y'\right] +n(n +1)y = 0 \]
Consider \(\,y'' +\lambda y =0~\) subject to the periodic boundary condition \(~y(-L)=y(L)\), \(\,y'(-L)=y'(L)\). \(\,\)Show that the eigenfunctions are
\[ \left\{ 1,\cos\frac{\pi}{L}x, \cos\frac{2\pi}{L}x, \cdots, \sin\frac{\pi}{L}x, \sin\frac{2\pi}{L}x, \,\cdots \right\}\]
This set, \(\,\)which is orthogonal on \(\,[-L,L]\), \(\,\)is the basis for the Fourier Series
\(~\)
The parametric Bessel differential equation is
\[ \frac{d}{dx}\left[xy'\right] +\left(\alpha^2x -\frac{n^2}{x}\right)y=0 \]
in which \(~r(x)=x\), \(~p(x)=x\), \(~q(x)=-n^2/x\), \(\,\)and \(~\lambda=\alpha^2\). \(\,\)The general solution of this equation is
\[ \color{blue}{y=c_1J_n(\alpha x)} +c_2Y_n(\alpha x) \]
Now \(r(0)=0\), \(\,\)and of two solutions \(J_n(\alpha x)\) and \(Y_n(\alpha x)\), \(\,\) only \(J_n(\alpha x)\) is bounded at \(x=0\)
The eigenvalues \(\lambda_i=\alpha_i^2, \;i=1,2,3,\cdots,\) \(\,\)are defined by means of a boundary condition at \(x=b\):
\[A_2J_n(\alpha b) +B_2\alpha J_n'(\alpha b) = 0 \]
For any choice of \(A_2\) and \(B_2\), \(\,\)not both zero, \(\,\)it is known that the above boundary condition gives an infinite number of roots \(x_i=\alpha_i b\). The eigenvalues are then \(\lambda_i=\alpha_i^2=(x_i/b)^2\)
Differential Recurrence Relations
\[ \begin{aligned} xJ_n'(x)&=xJ_{n-1}(x) -nJ_n(x)\\ xJ_n'(x)&=nJ_n(x)-xJ_{n+1}(x)\\ &\Downarrow \\ \frac{d}{dx}\left[ x^n J_n(x) \right ]&= x^n J_{n-1}(x)\\ \frac{d}{dx}\left[ x^{-n} J_n(x) \right ]&=-x^{-n} J_{n+1}(x) \end{aligned}\]
Square Norm
\[ \begin{aligned} \frac{d}{dx}\left[ xy' \right] &+\left( \alpha^2x -\frac{n^2}{x} \right)y = 0\\ &\Downarrow \;\times \;2xy' \\ \frac{d}{dx}\left[ xy' \right]^2 &+\left( \alpha^2x^2 -n^2 \right) \frac{d}{dx}[y]^2 =0\\ &\Downarrow {\tiny \text{integrating by parts on}\; [0,b]}\\ \left[ \left[xy' \right]^2 +(\alpha^2x^2 -n^2)y^2 \right]_0^b &=2\alpha^2 \int_0^b xy^2\,dx \\ &\Downarrow \; {\tiny y=J_n(\alpha x),\;y'=\alpha J_n'(\alpha x), \;J_n(0)=0\,\;\text{for}\; n>0} \\ \color{blue}{2\alpha^2 \int_0^b xJ_n^2(\alpha x)\,dx}\; &\color{blue}{= \alpha^2b^2 \left[ J_n'(\alpha b) \right]^2 +\left(\alpha^2b^2 -n^2\right)\left[ J_n(\alpha b) \right]^2} \end{aligned}\]
Consider three cases of the boundary condition
\[A_2J_n(\alpha b) +B_2\alpha J_n'(\alpha b) = 0\]
\(~\)
Case I : \(~~J_n(\alpha b)=0\;\) \(\Rightarrow\) \(\;\lambda_i=\alpha_i^2 =(x_i/b)^2 >0\)
\[ \begin{aligned} 2\alpha_i^2 \int_0^b xJ_n^2(\alpha_i x)\,dx &= \alpha_i^2b^2 \left[ J_n'(\alpha_i b) \right]^2 +(\alpha_i^2b^2 -n^2)\left[ J_n(\alpha_i b) \right]^2 \\ &\Downarrow \;xJ_n'(x)=n\underbrace{J_n(x)}_{=0} -xJ_{n+1}(x) \\ \color{blue}{\| J_n(\alpha_i x) \|^2 = \int_0^b xJ_n^2(\alpha_i x)}&\color{blue}{\,dx\; =\frac{b^2}{2} J_{n+1}^2(\alpha_i b)} \end{aligned}\]
\(~\)
from scipy.special import jv, jn_zeros
from mpl_toolkits.axisartist.axislines import SubplotZero
def make_plot_form():
fig = plt.figure(figsize=(7, 8))
ax1 = SubplotZero(fig, 211)
ax2 = SubplotZero(fig, 212)
fig.add_subplot(ax1)
ax1.axis["left", "right", "bottom", "top"].set_visible(False)
ax1.axis["xzero", "yzero"].set_visible(True)
ax1.axis["xzero", "yzero"].set_axisline_style("-|>")
ax1.set_xticks([2, 4, 6, 8, 10])
ax1.text(12, -0.03, r'$\alpha$', fontsize=14)
fig.add_subplot(ax2)
ax2.axis["left", "right", "bottom", "top"].set_visible(False)
ax2.axis["xzero", "yzero"].set_visible(True)
ax2.axis["xzero", "yzero"].set_axisline_style("-|>")
ax2.set_xticks([0.2, 0.4, 0.6, 0.8, 1.0])
ax2.text(1.1, -0.03, r'$x$', fontsize=14)
r1 = np.linspace(0, 11, 100)
r2 = np.linspace(0, 1, 100)
return ax1, ax2, r1, r2
def make_plot(case_no, ax1, ax2, r1, r2, mm, bc_txt, eigfn_txt):
alpha_0m = eval_alpha(mm)
y1 = eval_bc(r1)
ax1.plot(r1, y1, label=bc_txt)
ax1.plot(alpha_0m, np.zeros(mm), 'ro')
for m in range(mm):
ax1.text(alpha_0m[m] -0.2, 0.1, rf'$\alpha_{m +1}$',
fontsize=12)
ax1.legend(bbox_to_anchor=(1, 1), loc="upper left")
x_zeros = eval_x_zeros(mm)
for m in range(mm):
y2 = eval_eigfn(alpha_0m[m]*r2)
ax2.plot(r2, y2, label=eigfn_txt +rf'$(\alpha_{m +1}x)$')
if case_no == 1:
ax2.plot(x_zeros[:m +1] /alpha_0m[m], np.zeros(m +1), 'x')
elif case_no == 2:
ax2.plot(x_zeros[:m] /alpha_0m[m], np.zeros(m), 'x')
elif case_no == 3:
if m != 0:
ax2.plot(x_zeros[:m] /alpha_0m[m], np.zeros(m), 'x')
ax2.legend(bbox_to_anchor=(1, 1), loc="upper left")def eval_eigfn(x):
return jv(0, x)
def eval_bc(alpha):
return eval_eigfn(alpha)
def eval_alpha(mm):
return jn_zeros(0, mm)
def eval_x_zeros(mm):
return eval_alpha(mm)
def example_plot4():
ax1, ax2, r1, r2 = make_plot_form()
ax1.set_ylim(-0.5, 1.1)
ax1.set_yticks([-0.5, 0.5, 1.0])
ax2.set_ylim(-0.5, 1.1)
ax2.set_yticks([-0.5, 0.5, 1.0])
mm = 3
make_plot(1, ax1, ax2, r1, r2, mm, r'$J_0(\alpha)$', '$J_0$')
plt.show()
example_plot4()
def eval_eigfn(x):
return jv(1, x)
def eval_bc(alpha):
return eval_eigfn(alpha)
def eval_alpha(mm):
return jn_zeros(1, mm)
def eval_x_zeros(mm):
return eval_alpha(mm)
def example_plot4_():
ax1, ax2, r1, r2 = make_plot_form()
ax1.set_ylim(-0.5, 1.1)
ax1.set_yticks([-0.5, 0.5, 1.0])
ax2.set_ylim(-0.5, 1.1)
ax2.set_yticks([-0.5, 0.5, 1.0])
mm = 3
make_plot(1, ax1, ax2, r1, r2, mm, r'$J_1(\alpha)$', '$J_1$')
plt.show()
example_plot4_()
\(~\)
Case II : \(\,\)(\(h\geq 0\))
\[~ hJ_n(\alpha b) +\alpha bJ_n'(\alpha b)=0 \;\Rightarrow \;\lambda_i=\alpha_i^2 =\left(\frac{x_i}{b}\right)^2 \;\text{ except for }\; h=0\; \text{ and } \; n=0 \]
\[ \begin{aligned} 2\alpha_i^2 \int_0^b xJ_n^2(\alpha_i x)\,dx &= \alpha_i^2b^2 \left[ J_n'(\alpha_i b) \right]^2 +(\alpha_i^2b^2 -n^2)\left[ J_n(\alpha_i b) \right]^2 \\ \\ &\Downarrow \;\alpha_ibJ_n'(\alpha_i b)=-hJ_n(\alpha_i b) \\ \\ \color{blue}{\| J_n(\alpha_i x) \|^2 = \int_0^b xJ_n^2(\alpha_i x)\,dx} \; & \color{blue}{=\frac{\alpha_i^2b^2 -n^2 +h^2}{2\alpha_i^2} J_n^2(\alpha_i b) } \end{aligned}\]
\(~\)
from scipy.special import j0, j1
from scipy.optimize import fsolve
def eval_eigfn(x):
return j0(x)
def eval_bc(alpha):
return j0(alpha) -alpha *j1(alpha)
def eval_alpha(mm):
alpha_0m = np.zeros(mm)
for i in range(mm):
alpha_0m[i] = fsolve(eval_bc, 1 +3*i)
return alpha_0m
def eval_x_zeros(mm):
return jn_zeros(0, mm -1)
def example_plot4__():
ax1, ax2, r1, r2 = make_plot_form()
ax1.set_ylim(-4, 4)
ax1.set_yticks([-3, -2, -1, 1, 2, 3])
ax2.set_ylim(-0.5, 1.2)
ax2.set_yticks([-0.5, 0.5, 1.0])
mm = 4
make_plot(2, ax1, ax2, r1, r2, mm,
r"$J_0(\alpha) +\alpha J_0'(\alpha)$", '$J_0$')
plt.show()
example_plot4__()
\(~\)
Case III : \(\,\)(\(h=0\) and \(n=0\)) : \(~~J_0'(\alpha b)=0\;\)
It is the only situation for which \(\lambda=0\) is an eigenvalue
\[ \begin{aligned} \color{blue}{\alpha_1 = 0 \rightarrow \| J_0(\alpha_1 x) \|^2}\; &\color{blue}{= \int_0^b x\,dx =\frac{b^2}{2}} \\ \\ 2\alpha_i^2 \int_0^b xJ_0^2(\alpha_i x)\,dx &= \alpha_i^2b^2 \left[ J_0'(\alpha_i b) \right]^2 +(\alpha_i^2b^2 -0^2)\left[ J_0(\alpha_i b) \right]^2 \\ &\Downarrow \\ \color{blue}{\alpha_i,\; i=2,3,4,\cdots \rightarrow \| J_0(\alpha_i x) \|^2} \; &\color{blue}{ =\int_0^b xJ_0^2(\alpha_i x)\,dx =\frac{b^2}{2}J_0^2(\alpha_i b)} \end{aligned}\]
\(~\)
from scipy.special import jnp_zeros
def eval_eigfn(x):
return j0(x)
def eval_bc(alpha):
return -j1(alpha)
def eval_alpha(mm):
return np.append([0], jnp_zeros(0, mm -1))
def eval_x_zeros(mm):
return jn_zeros(0, mm -1)
def example_plot4___():
ax1, ax2, r1, r2 = make_plot_form()
ax1.set_ylim(-0.8, 0.5)
ax1.set_yticks([-0.5, 0.5])
ax2.set_ylim(-0.5, 1.2)
ax2.set_yticks([-0.5, 0.5, 1.0])
mm = 4
make_plot(3, ax1, ax2, r1, r2, mm, r"$J_0'(\alpha)$", '$J_0$')
plt.show()
example_plot4___()
\(~\)
The set \(\{J_n(\alpha_i x), \;i=0,1,2,\cdots\}\) is orthogonal with respect to the weight function \(p(x)=x~\) on \(~[0, b]\)
\[\int_{0}^b x\, J_n(\alpha_i x) \, J_n(\alpha_j x)\,dx=0,\; i\neq j\]
The Fourier-Bessel series of a function \(\,f\) defined on the interval \((0, b)\) is given by
\[ f(x)=\sum_{i=1}^\infty c_i J_n(\alpha_i x)\]
where
\[ c_i=\frac{\displaystyle\int_0^b x\,f(x)\,J_n(\alpha_i x)\,dx}{\displaystyle\int_0^b x \, J_n^2(\alpha_i x)\,dx} =\color{blue}{\frac{1}{\| J_n(\alpha_i x)\|^2}}\int_0^b x\,f(x) \, J_n(\alpha_i x)\,dx \]
\(~\)
Legendre’s differential equation is
\[(1-x^2)y'' -2xy'+n(n+1)y=0\;\]
\[\text{or}\]
\[\frac{d}{dx}\left[ (1-x^2)y' \right]+n(n+1)y=0\]
in which \(~r(x)=1-x^2, \,p(x)=1, \,\lambda=n(n+1),\) and \(~q(x)=0\)
The Legendre polymomials \(P_n(x)\) are the only solutions that are bounded on the interval \([\text{-}1, \,1]\)
\(~\)
from scipy.special import eval_legendre
fig = plt.figure(figsize=(5, 5))
ax = SubplotZero(fig, 111)
fig.add_subplot(ax)
ax.axis["left", "right", "bottom", "top"].set_visible(False)
ax.axis["xzero", "yzero"].set_visible(True)
ax.axis["xzero", "yzero"].set_axisline_style("-|>")
ax.set_xticks([-1, -0.5, 0.5, 1], ['-1', '', '', '1'])
ax.set_yticks([-1, -0.5, 0.5, 1], ['-1', '', '', '1'])
ax.grid(ls=':')
ax.text(-0.1, 1.3, r'$P_n(x)$', fontsize=14)
ax.text(1.3, -0.05, r'$x$', fontsize=14)
ax.set_xlim(-1.2, 1.2)
ax.set_ylim(-1.2, 1.2)
x = np.linspace(-1, 1, 100)
for n in range(5):
y = eval_legendre(n, x)
ax.plot(x, y, label=f'$P_{n}(x)$')
ax.legend(loc='lower left', bbox_to_anchor=(1, 0))
plt.show()
\(~\)
The set \(\{P_n(x), \;n=0,1,2,\cdots\}\) is orthogonal with respect to the weight function \(p(x)=1\) on \([\text{-}1,\,1]\)
\[\int_{-1}^1 P_m(x)P_n(x)\,dx=0,\; m\neq n\]
The Fourier-Legendre series of a function \(\,f\) defined on the interval \((\text{-}1,1)\) is given by
\[ f(x)=\sum_{n=0}^\infty c_n P_n(x)\]
where
\[ c_n=\frac{\displaystyle\int_{-1}^1 f(x)P_n(x)\,dx}{\displaystyle\int_{-1}^1 P_n^2(x)\,dx} =\color{blue}{\frac{2n+1}{2}}\int_{-1}^1 f(x)P_n(x)\,dx \]
\(~\)
Hermite
\[y'' -2x y' +2n y=0\]
Laguerre
\[xy'' +(1-x)y' +ny=0\]
Chebyshev
\[(1-x^2) y'' -xy' +n^2y = 0\]
Hypergeometric
\[x(1-x)y'' +[c - (a + b +1)x y'] -ab y = 0\]
\(~\)