Examples on how to integrate a computational perspective in physics courses

1. Ordinary differential equations (ODE): RLC circuit
2. ODE: Classical pendulum
3. ODE: Solar system
4. and many more cases

Classical pendulum with damping and external force as it could appear in a mechanics course (PHY 321) $$ml\frac{d^2\theta}{dt^2}+\nu\frac{d\theta}{dt} +mgsin(\theta)=Acos(\omega t).$$ Easy to solve numerically and then visualize the solution. Almost the same equation for an RLC circuit in the electromagnetism course (PHY 481/482) $$L\frac{d^2Q}{dt^2}+\frac{Q}{C}+R\frac{dQ}{dt}=Acos(\omega t).$$

Classical pendulum equations with damping and external force $$\frac{d\theta}{d\hat{t}} =\hat{v},$$ and $$\frac{d\hat{v}}{d\hat{t}} =Acos(\hat{\omega} \hat{t})-\hat{v}\xi-\sin(\theta),$$ with $\omega_0=\sqrt{g/l}$, $\hat{t}=\omega_0 t$ and $\xi = mg/\omega_0\nu$.

The RLC circuit $$\frac{dQ}{d\hat{t}} =\hat{I},$$ and $$\frac{d\hat{I}}{d\hat{t}} =Acos(\hat{\omega} \hat{t})-\hat{I}\xi-Q,$$ with $\omega_0=1/\sqrt{LC}$, $\hat{t}=\omega_0 t$ and $\xi = CR\omega_0$.

The equations are essentially the same. Great potential for abstraction. Use the same program with small modifications.

These physics examples can all be studied using almost the same types of algorithms, simple eigenvalue solvers or Gaussian elimination with almost the same starting matrix!

1. A buckling beam and tridiagonal Toeplitz matrices (mechanics and mathematical methods), eigenvalue problems
2. A particle in an infinite potential well, quantum eigenvalue problems
3. A particle (or two) in a general quantum well, quantum eigenvalue problems
4. Poisson's equation in one dim, linear algebra (electromagnetism)
5. The diffusion equation in one dimension (Statistical Physics), linear algebra
6. and many other cases

This is a two-point boundary value problem $$R \frac{d^2 u(x)}{dx^2} = -F u(x),$$ where $u(x)$ is the vertical displacement, $R$ is a material specific constant, $F$ the force and $x \in [0,L]$ with $u(0)=u(L)=0$.

Scale equations with $x = \rho L$ and $\rho \in [0,1]$ and get (note that we change from $u(x)$ to $v(\rho)$) $$\frac{d^2 v(\rho)}{dx^2} +K v(\rho)=0,$$ a standard eigenvalue problem with $K= FL^2/R$.

If you replace $R=-\hbar^2/2m$ and $-F=\lambda$, we have the quantum mechanical variant for a particle moving in a well with infinite walls at the endpoints.

Discretize the second derivative and the rhs $$-\frac{v_{i+1} -2v_i +v_{i-i}}{h^2}=\lambda v_i,$$ with $i=1,2,\dots, n$. We need to add to this system the two boundary conditions $v(0) =v_0$ and $v(1) = v_{n+1}$. The so-called Toeplitz matrix (special case from the discretized second derivative) $$\mathbf{A} = \frac{1}{h^2}\begin{bmatrix} 2 & -1 & & & & \\ -1 & 2 & -1 & & & \\ & -1 & 2 & -1 & & \\ & \dots & \dots &\dots &\dots & \dots \\ & & &-1 &2& -1 \\ & & & &-1 & 2 \\ \end{bmatrix}$$ with the corresponding vectors $\mathbf{v} = (v_1, v_2, \dots,v_n)^T$ allows us to rewrite the differential equation including the boundary conditions as a standard eigenvalue problem $$\mathbf{A}\mathbf{v} = \lambda\mathbf{v}.$$ The Toeplitz matrix has analytical eigenpairs!! Adding a potential along the diagonals allows us to reuse this problem for many types of physics cases.

Assume we want to solve the radial part of Schroedinger's equation for one particle in three dimensions. This equation reads $$-\frac{\hbar^2}{2 m} \left ( \frac{1}{r^2} \frac{d}{dr} r^2 \frac{d}{dr} - \frac{l (l + 1)}{r^2} \right )R(r) + V(r) R(r) = E R(r).$$ Suppose in our case $V(r)$ is the harmonic oscillator potential $(1/2)kr^2$ with $k=m\omega^2$ and $E$ is the energy of the harmonic oscillator in three dimensions.

Scale now the equations with $\rho = r/\alpha$ where $\alpha$ is a constant of dimension length.

Manipulating the equations by requiring $$\frac{mk}{\hbar^2} \alpha^4 = 1,$$ which define defines a natural length scale (like the Bohr radius does) $$\alpha = \left(\frac{\hbar^2}{mk}\right)^{1/4}.$$ Defining $$\lambda = \frac{2m\alpha^2}{\hbar^2}E,$$ we can rewrite Schroedinger's equation as $$-\frac{d^2}{d\rho^2} v(\rho) + \rho^2v(\rho) = \lambda v(\rho) .$$ This is similar to the equation for a buckling beam except for the potential term.

The last example shows the potential of combining numerical algorithms with analytical results (or eventually symbolic calculations). A simple change of potential gives a new physics case, example of a box potential

# Different types of potentials
def potential(r):
    if r >= 0.0 and r <= 10.0:
        V = -0.05
    else:
        V =0.0
    return V

This allows students and teachers to

• make abstraction and explore other physics cases easily where no analytical solutions are known
• Validate and verify their algorithms.
• Including concepts like unit testing, one has the possibility to test and validate several or all parts of the code.
• Validation and verification are then included naturally and one can develop a better attitude to what is meant with an ethically sound scientific approach.
• The above example allows the student to also test the mathematical error of the algorithm for the eigenvalue solver by changing the number of integration points. The students get trained from day one to think error analysis.
• The algorithm can be tailored to any kind of one-particle problem used in quantum mechanics or eigenvalue problems
• A simple rewrite allows for reuse in linear algebra problems for solution of say Poisson's equation in electromagnetism, or the diffusion equation in one dimension.
• With an ipython notebook the students can keep exploring similar examples and turn them in as their own notebooks.