Week 10 Day 3: ODE problems

Objectives

Discuss debuggers
Solve a BVP
Look at some other ODE problems

Debugging

Using a debugger has several parts:

Starting the debugger
- Classic way: import pdb; pdb.set_trace() (This is the standard library debugger)
- Python 3.7 way: breakpoint() (This will start the current debugger)
- IPython offers %debug to jump in after an error is raised! (like pdb.pm())
- You can also start at the beginning (see docs)
Controlling the debugger
- pyb (and ipydb) are sort of based on classic debugging tools - the syntax is a bit odd, but will be useful in other languages.
- A common term is the "stack", short in this case for the "call-stack", the series of functions that call each other to get to the current point.
  - Note that the "stack" is also a term for memory locations in compiled programs - and is unrelated to this usage.

short	name	action	Names in other debuggers
`h`	`help`	Print out help
`p`	`print`	Print out the value of a variable or expression
`u`	`up`	Go up one in the stack
`d`	`down`	Go down one in the stack
`w`	`where`	Show you where you are	Stack trace
`s`	`step`	Move forward one computation	Step in
`n`	`next`	More forward one line	Step over
`c`	`continue`	Continue on till next stop
`b`	`breakpoint`	Set a breakpoint somehere
`q`	`quit`	Quit

## A program to debug
def simple_function(div):
    simple_x = 2
    simple_y = div
    
    # import pdb; pdb.set_trace()
    
    return fancy_function(simple_x, simple_y)

def fancy_function(x, y):
    r = x / y
    return r

# simple_function(0)

# %debug

# simple_function(2)

# import pdb

# pdb.run('simple_function(0)')

Other tools: Tracing

Some IDEs, like PyCharm, may offer enhanced and partially graphical debuggers. And, in Python 3.7, the new built-in breakpoint() can call the fancy debugger instead of the builtin one!

Let's look at something similar: Tracing. If you have a complex piece of Python code and you want to have an idea of what the control flow looks like, you can trace it.

%%writefile simple.py

## A program to debug
def simple_function(div):
    simple_x = 2
    simple_y = div
    
    return fancy_function(simple_x, simple_y)

def fancy_function(x, y):
    r = x / y
    return r

def main():
    print(simple_function(2))

if __name__ == '__main__':
    main()

!python simple.py

!python -m trace --trace simple.py

!python -m trace --listfuncs simple.py

Boundary value problems (BPV)

Let's look at a square well:

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate

$$ \begin{align} \psi'' + \left( \frac{2 m}{\hbar^2} V_0 - \kappa^2 \right) \psi = 0,& \quad \textrm{for }|x|\le a \\ \psi'' - \kappa^2 \psi = 0,& \quad \textrm{for }|x| > a \end{align} $$

Again, we write this as a set of first order equations:

$$ u = \left(\begin{matrix}\psi' \\ \psi\end{matrix}\right) $$$$ u' = \left(\begin{matrix}\psi'' \\ \psi'\end{matrix}\right) = \left(\begin{matrix}\left( \frac{2 m}{\hbar^2} V_0 - \kappa^2 \right) \psi \\ \psi'\end{matrix}\right) $$

where we set $V_0$ to $0$ when $|x|$ is larger than $a$. Boundary conditions:

$$ \psi(x_\mathrm{max}) = e^{-\kappa x_\mathrm{max}} \\ \psi(-x_\mathrm{max}) = e^{-\kappa x_\mathrm{max}} $$

Double click to see broken code for solve_bvp (in a live notebook, not a viewer).

""" From "COMPUTATIONAL PHYSICS" & "COMPUTER PROBLEMS in PHYSICS"
    by RH Landau, MJ Paez, and CC Bordeianu (deceased)
    Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, 
    C Bordeianu, Univ Bucharest, 2017. 
    Please respect copyright & acknowledge our work."""

# QuantumEigen.py:             Finds E and psi via rk4 + bisection
 
# mass/((hbar*c)**2)= 940MeV/(197.33MeV-fm)**2 =0.4829
# well width=20.0 fm, depth 10 MeV, Wave function not normalized.  

import numpy as np
import matplotlib.pyplot as plt

eps       = 1e-3                                     # Precision
n_steps   = 501 
E         = -17                                       # E guess
h         = 0.04
count_max = 100
Emax      = 1.1*E                                     # E limits
Emin      = E/1.1            

def f(x, y, E):
    return np.array([y[1], -0.4829*(E-V(x))*y[0]])

def V(x):
    # Well depth
    return -16 if abs(x) < 10 else 0
    
def rk4(t, y, h, Neqs, E):
    ydumb = np.zeros((Neqs),float)
    k1 = np.zeros((Neqs),float)
    k2 = np.zeros((Neqs),float)
    k3 = np.zeros((Neqs),float)
    k4 = np.zeros((Neqs),float)

    F = f(t, y, E)
    k1 = h*F
    ydumb = y + k1/2
        
    F = f(t + h/2, ydumb, E)
    k2 = h*F
    ydumb = y + k2/2
    
    F = f(t + h/2, ydumb, E)
    k3 =  h*F
    ydumb = y + k3
        
    F = f(t + h, ydumb, E)
    k4 = h*F
    
    y += (k1 + 2*(k2 + k3) + k4)/6
        

def diff(E, h):
    i_match = n_steps//3                     # Matching radius
    nL = i_match + 1  
    
    y = np.array([1e-15, 1e-15*np.sqrt(-E*0.4829)]) # Initial left wf  
    
    for ix in range(0,nL + 1):
        x = h * (ix  -n_steps/2)
        rk4(x, y, h, 2, E)
        
    left = y[1]/y[0]                        # Log  derivative
    
    y[0] = 1e-15;                    # For even;  reverse for odd
    y[1] = -y[0]*np.sqrt(-E*0.4829)   # Initialize R wf
    
    for ix in range(n_steps, nL+1, -1):
        x = h*(ix + 1 - n_steps/2)
        rk4(x, y, -h, 2, E)
        
    right = y[1]/y[0]                 # Log derivative
    return (left - right) / (left + right) 


def plot(E, h):
    # Repeat integrations for plot
    x = 0. 
    n_steps = 1501                       # integration steps
    
    y = np.zeros((2),float)
    yL = np.zeros((2,505),float) 
    i_match = 500                             # Matching point
    nL = i_match + 1
    
    y[0] = 1e-40                            # Initial left wf
    y[1] = -np.sqrt(-E*0.4829) * y[0]
    
    for ix in range(nL+1):                          
        yL[0][ix] = y[0]
        yL[1][ix] = y[1]
        x = h * (ix -n_steps/2)
        rk4(x, y, h, 2, E)
    
    y[0] = -1.E-15                # For even;  reverse for odd
    y[1] = -np.sqrt(-E*0.4829)*y[0]
    j = 0
    
    Rwf_x = np.zeros(n_steps - 3 - nL)
    Rwf_y = np.zeros(n_steps - 3 - nL)
    for ix in range(n_steps-1, nL+2, -1):          # right WF
        x = h * (ix + 1 -n_steps/2)                # Integrate in
        rk4(x, y, -h, 2, E)
        Rwf_x[j] = 2*(ix + 1 - n_steps/2)
        Rwf_y[j] = y[0]*35e-9
        j += 1
    x = x-h              
    normL = y[0]/yL[0][nL]
    j=0
    
    # Renormalize L wf & derivative
    Lwf_x = np.zeros(nL+1)
    Lwf_y = np.zeros(nL+1)
    for ix in range(nL+1):
        x = h * (ix-n_steps/2 + 1) 
        y[0] = yL[0][ix]*normL 
        y[1] = yL[1][ix]*normL
        Lwf_x[j] = 2*(ix - n_steps/2 + 1)
        Lwf_y[j] = y[0]*35e-9           # Factor for scale
        j +=1
        
    print(E)
    plt.plot(Rwf_x, Rwf_y)
    plt.plot(Lwf_x, Lwf_y)
    plt.show()
        
for count in range(0,count_max+1):
    # Iteration loop
    E = (Emax + Emin)/2                     # Divide E range
    Diff = diff(E, h) 
    
    if (diff(Emax, h)*Diff > 0):
        Emax = E  # Bisection algor
    else:
        Emin = E
        
    if abs(Diff) < eps:
        break
    
    plot(E, h)
    
        
print("Final eigenvalue E = ",E)
print("iterations, max = ",count)

Problem 2: Projectile motion with drag

We can select a power $n=1$ or $2$. The friction coefficient is $k$.

$$ \ddot{x} = - k \dot{x}^n \frac{\dot{x}}{v} \\ \ddot{y} = - g - k \dot{y}^n \frac{\dot{y}}{v} \\ v = \sqrt{\dot{x}^2 + \dot{y}^2} $$

We can rewrite it:

$$ u = \left(\begin{matrix} x \\ \dot{x} \\ y \\ \dot{y} \end{matrix}\right) $$

Then:

$$ \dot{u} = \left(\begin{matrix} \dot{x} \\ - k \dot{x}^n \frac{\dot{x}}{v} \\ \dot{y} \\ - g - k \dot{y}^n \frac{\dot{y}}{v} \end{matrix}\right) = \left(\begin{matrix} u_1 \\ - k u_1^n \frac{u_1}{\sqrt{u_1^2 + u_3^2}} \\ u_3 \\ - g - k u_3^n \frac{u_3}{\sqrt{u_1^2 + u_3^2}} \end{matrix}\right) $$

We can set the initial positions to 0 and give it a starting velocity, and compare no friction to friction versions.

import scipy.integrate

v0 = 22
angle = 34
g = 9.8
k = 0.8
n = 1

v0x = v0*np.cos(np.radians(angle))
v0y = v0*np.sin(np.radians(angle))

t_eval = np.linspace(0,2)
analytic_x = v0x*t_eval
analytic_y = v0y*t_eval - 0.5*g*t_eval**2

def f(t, u):
    v = np.sqrt(u[1]**2 + u[3]**2)
    return np.stack([
        u[1],
        -k * u[1]**n * u[1] / v,
        u[3],
        -g - k * u[3]**n * u[3] / v
    ])

res = scipy.integrate.solve_ivp(f, [0, 2], [0, v0x, 0, v0y], t_eval=t_eval)
print(res.message)

plt.figure(figsize=(8,3))
plt.plot(res.y[0], res.y[2])
plt.plot(analytic_x, analytic_y)

""" From "COMPUTATIONAL PHYSICS" & "COMPUTER PROBLEMS in PHYSICS"
    by RH Landau, MJ Paez, and CC Bordeianu (deceased)
    Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, 
    C Bordeianu, Univ Bucharest, 2017. 
    Please respect copyright & acknowledge our work."""

# ProjectileAir.py: Order dt^2 projectile trajectory + drag

import numpy as np
import matplotlib.pyplot as plt

v0 = 22
angle = 34.
g = 9.8
kf = 0.8
N = 20

v0x = v0*np.cos(angle*np.pi/180)
v0y = v0*np.sin(angle*np.pi/180)
T = 2*v0y/g
H = v0y**2/2/g
R = 2*v0x*v0y/g

print('No Drag T =', T, ', H =', H, ', R =', R)

def plotNumeric(k):
    vx = v0*np.cos(angle*np.pi/180)
    vy = v0*np.sin(angle*np.pi/180)
    x = 0.0
    y = 0.0
    dt =  vy/g/N*1.5
    
    print("\n       With Friction  ")
    print("       x            y")
    
    xy = np.empty((2, N))
    xy[:,0] = 0,0

    for i in range(N-1):
        vx = vx - k*vx*dt
        vy = vy - g*dt - k*vy*dt
        x = x + vx*dt
        y = y + vy*dt
        xy[:,i+1] = x, y
        print(" %13.10f  %13.10f " % (x,y))
    
    plt.plot(*xy)

def plotAnalytic():
    v0x = v0*np.cos(angle*np.pi/180)
    v0y = v0*np.sin(angle*np.pi/180)
    dt = v0y/g/N*1.8
    print("\n       No Friction  ")
    print("        x            y")
    
    xy = np.empty((2, N))
    for i in range(N):
        t = i*dt
        x = v0x*t
        y = v0y*t - g*t*t/2
        xy[:,i] = x, y
        print(" %13.10f  %13.10f" % (x ,y))
        
    plt.plot(*xy)
        
plotNumeric(kf)
plotAnalytic()