How do you even teach series solutions without Jupyter?
Time to talk Bessel!
First, let's set up our basic environment
import numpy as np
from matplotlib import pyplot as plt
from scipy.special import gamma
%matplotlib inline
Now, let's define the solutions to Bessel's equation. We have (12.9)
$J_p(x) = \sum_{n=0}^\infty \frac{(-1)^n}{\Gamma(n+1)\Gamma(n+1+p)} \left(\frac{x}{2} \right)^{2n+p}$
And we'll do our normal thing, whwere we make up a bunch of points upon which to calculate our function, then start with zero and add terms.
First, a bit of programming syntax, +=
thing = 2
thing = thing + 2
print(thing)
thingything = 6
thingything += 2 # += means thing = thing + 2.
print(thingything)
Now, let's define a function to calculate bessel functions, just exactly implementing the mathematical definition above
def Jp(p,x,N):
result = 0
for n in range(0,N+1):
#result = result + new stuff, for term n
result += ((x/2)**(2*n+p)) * ((-1)**n)/(gamma(n+1)*gamma(n+1+p))
return result
Jp(4,3,4)
Cool. That gives us back the value at one point. Using numpy arrays, we can get a range of values, so that we can plot them.
def Jp(p,x,N):
result = np.zeros_like(x)
for n in range(0,N+1):
#result = result + new stuff, for term n
result += ((x/2)**(2*n+p)) * ((-1)**n)/(gamma(n+1)*gamma(n+1+p))
return result
x = np.linspace(0,10,1000)
plt.plot(x,Jp(4,x,7))
Make it interactive¶
Here's one of the things that's just amazingly useful about Jupyter notebooks. I'd like to interact with teh above thing.
from ipywidgets import interact, fixed
def plotJ(p,N):
x = np.linspace(0,10,1000)
plt.plot(x,Jp(p,x,N))
plotJ(4,9)
interact(plotJ,p=(-5,5),N=(0,30))
Oh, now I'd like to mess around with the x axis, so I can see what happens farther out.
def plotJ(p=4,N=10,xmax=10):
x = np.linspace(0,xmax,1000)
plt.plot(x,Jp(p,x,N))
interact(plotJ,p=(-5,5),N=(0,300,10),xmax=(10,100,10))
We learned a few interesting lessons from the above, especially about numerical issues. Among them
Adding more terms doesn't always help. It sometimes leads to numerical instability
We also noticed that sometimes the xmax parameter didn't work. Here's an example, where we should be going from 0 to 60, but it looks like we only go from 0 to about 37:
xmax=60
x = np.linspace(0,xmax,1000)
jpx = Jp(p=4,x=x,N=120)
plt.plot(x,jpx)
Our first hint is that we're getting a bunch of errors about division. Let's look at the data, but maybe with only 100 data points:
xmax=60
x = np.linspace(0,xmax,100)
jpx = Jp(p=4,x=x,N=120)
plt.plot(x,jpx)
jpx
All of those "nan"s at the end are what happened when we tried to divide by zero. By default, those get ignored during plotting. Here, ignoring them means stopping the plot early. Cool.