Exercises on differential equation

The pendulum equation

We consider the equation

\theta^{\prime\prime} = -\sin{\theta}

Write a function that returns an array containing the position and the angular velocity of the pendulum for N instants t_i between 0 and T.

The initial position is \theta=0. Plot the phase space trajectory for different values of the initial velocity. Angle will be represented between -π and π.

Solution

from scipy.integrate import ode
from numpy import *

def f(t, Y):
    y, yprime=Y
    return array([yprime, -sin(y)])

def pendule(theta_0, vit_ang, T, N=1000):
    r = ode(f).set_integrator('vode')
    r.set_initial_value(array([theta_0, vit_ang]),0)
    TT = linspace(0,T,N+1)
    output = [[0, theta_0, vit_ang]] # Output is a list of list
    for i, t in enumerate(TT[1:]):
        r.integrate(t)
        output.append([t, r.y[0], r.y[1]])
    return array(output) # This is a 2D array

figure(1)
clf()
Tvit_ini = linspace(-3,3, 61)
for vit_ini in Tvit_ini:
    a = pendule(0, vit_ini, 100)
    a[:,1] = ((a[:,1]+pi)%(2*pi) - pi) # Tricks to be between -pi and pi
    plot(a[:,1], a[:,2], '+')

xlim(-pi,pi)

Solving the Schrödinger equation using the finite element method

Let us consider the Schrödinger equation with hbar = m = 1

\frac{d\psi}{dt} =-i\left( -\frac12 \frac{d^2\psi}{dx^2} + V(x)\psi\right)

The potential is V(x) = \frac12 k x^2 with \kappa=.5.

To solve the equation, we will truncate the x-axis to values between x_\mathrm{min} and x_\mathrm{max}. We will also discretize the x-axis with small steps (\Delta x).

The term \frac{d^2\psi}{dx^2} will be approximated using \frac{\psi(x+\Delta x) - 2\psi(x) + \psi(x-\Delta x)}{\Delta x^2}.

For the initial state, we will take a Gaussian distribution e^{-\alpha(x-x_0)^2}. We will use \alpha = \frac 12 and x_0=1.

Calculate and plot |{\psi(x,t)}|^2 as a function of x for t=1 using the zvode solver with an absolute precision of 10{^-3}.

Solution

from scipy.integrate import ode
from numpy import *

Deltax = 0.1
X_MIN = -10
X_MAX = 10

alpha = .4
x0 = 1

kappa = .5
T = 10
N = 1000

x = linspace(X_MIN, X_MAX, (X_MAX-X_MIN)/Deltax + 1)

psi_0 = exp(-alpha*(x-x0)**2)
psi_0 = psi_0 / sqrt(sum(abs(psi_0**2)))

V = -.5*kappa*x**2

psi_dot = zeros(len(x), dtype="complex128")

def f(t, psi):
    energy = .5*(psi[2:] - 2*psi[1:-1] + psi[:-2])/Deltax**2
    pot = V[1:-1]*psi[1:-1]
    psi_dot[1:-1] = -1j*(energy + pot)
    return psi_dot

r = ode(f).set_integrator('zvode', atol=1E-3)
r.set_initial_value(psi_0,0)

psi_f = r.integrate(1)

figure(0)
plot(x, abs(psi_f)**2)

# Calculate the mean position and relative width.
figure(1)
r = ode(f).set_integrator('zvode', atol=1E-3)
r.set_initial_value(psi_0,0)

dT = linspace(0,T,N+1)
out = []
for i, t in enumerate(dT):
    if t>r.t:
        r.integrate(t)
    print i
    out.append([t, sum(x*abs(r.y)**2)/sum(abs(r.y)**2), (mean(x**2*abs(r.y)**2) - mean(x*abs(r.y)**2)**2)/mean(abs(r.y)**2)])

Table Of Contents

Previous topic

Exercises on fit

Next topic

Exercises on Fourier analysis

This Page