Python Exercise 2

Fredrik Eikeland Fossan

Apr 19, 2018

Newton's first differential equation

One of the earliest known differential equations, which Newton solved with series expansion in 1671 is: $$ \begin{align} y'(x) =1-3x+y+x^2+xy,\ y(0)=0 \label{eq:newton} \end{align} $$

Newton gave the following solution: $$ \begin{align} y(x) \approx x-x^2+\frac{x^3}{3}-\frac{x^4}{6}+ \frac{x^5}{30}-\frac{x^6}{45} \label{eq:newton_approx} \end{align} $$

Today it is possible to give the solution on closed form with known functions as follows, $$ \begin{align} y(x)=&3\sqrt{2\pi e}\cdot \exp\left[x\left(1+\frac{x}{2}\right)\right]\cdot \left[\text{erf}\left(\frac{\sqrt{2}}{2}(1+x)\right)-\text{erf}\left(\frac{\sqrt{2}}{2}\right)\right] + 4\cdot\left[1-\exp[x\left(1+\frac{x}{2}\right)\right]-x \label{eq:newton_analytical} \end{align} $$

Note the combination $ \sqrt{2\pi e} $.

Problem 1: Euler's method

Solve Eq. \eqref{eq:newton} using Euler's method. Plot and compare with Newton's solution Eq. \eqref{eq:newton_approx} and the analytical solution Eq. \eqref{eq:newton_analytical}. Plot between $ x=0 $ and $ x=1.5 $

Hint 1.

The function scipy.special.erf will be needed. See http://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.special.erf.html.

Hint 2.

Figure should look something like this:

Figure 1: Problem 1 with 101 samplepoints between $ x=0 $ and $ x=1.5 $ .

Solution.

# python_exercises/pexercise2/problem1.py

import numpy as np
import matplotlib.pylab as plt
from scipy.special import erf


def euler(func, y_0, time):
    """ Generic implementation of the euler scheme for solution of systems of ODEs (note that x may be interchanged by e.g. x):
            y0' = y1
            y1' = y2
                .
                .
            yN' = f(yN-1,..., y1, y0, t)
            
            method:
            y0^(n+1) = y0^(n) + h*y1
            y1^(n+1) = y1^(n) + h*y2
                .
                .
            yN^(n + 1) = yN^(n) + h*f(yN-1, .. y1, y0, t)

        Args:
            func(function): func(y, t) that returns y' at time t; [y1, y2,...,f(yn-1, .. y1, y0, t)] 
            y_0(array): initial conditions
            time(array): array containing the time to be solved for
        
        Returns:
            y(array): array/matrix containing solution for y0 -> yN for all timesteps"""

    y = np.zeros((np.size(time), np.size(y_0)))
    y[0,:] = y_0

    for i in range(len(time)-1):
        h = time[i+1] - time[i]
        y[i+1,:]=y[i,:] + np.asarray(func(y[i,:], time[i]))*h

    return y

def newton_func(y, x):
    """ function that returns the RHS of the Newton ODE:
    
        y' = 1 - 3x + y + x^2 + xy
        
        Args:
            y(array): array [y0, y1] at distance x
            x(float): current position
        
        Returns:
            dy(float): y' 
        """
        
    dy = 1 - 3*x + y + x**2 + x*y
    return dy

x = np.linspace(0, 1.5, 101) #allocate Space
y_0 = 0

Y = euler(newton_func, y_0, x)
Y_newton = x - x**2 + x**3/3 - x**4/6 + x**5/30 - x**6/45

e = np.exp(1)
pi = np.pi

Y_analytical = 3*np.sqrt(2*pi*e)*np.exp(x*(1 + x/2))*(erf(np.sqrt(2)*(1 + x)/2) - erf(np.sqrt(2)/2)) + 4*(1-np.exp(x*(1 + x/2))) - x

plt.figure()
plt.plot(x, Y_analytical)
plt.plot(x, Y)
plt.plot(x, Y_newton)
plt.legend(['analytical', 'euler', 'Newton'], loc='best', frameon=False)
plt.xlabel('x')
plt.ylabel('y')
#plt.savefig('problem1.png', transparent=True)
plt.show()

Mathematical pendulum with friction

Figure 2: Mathematical pendulum.

Figure 2 shows a mathematical pendulum. By modeling the friction forces acting on the pendulum as $ \mu \frac{d \theta}{dt} $, the motion can be described by the following ODE: $$ \begin{align} \frac{d ^2\theta}{dt^2} + \frac{\mu}{m} \frac{d \theta}{dt} + \frac{g}{l} \sin\theta &=0, \label{eq:1} \end{align} $$

where $ \mu $ is a friction coefficient, $ l $ is the length og the rod and $ m $ is the mass of the pendulum with initial conditions $$ \begin{align} &\theta(0) =\theta_0 \label{_auto1}\\ &\frac{d \theta}{dt}(0)=\dot\theta_0. \label{_auto2} \end{align} $$

Assume that the pendulum wire is a massless rod, such that $ -\pi \leq \theta_0 \leq \pi $. Unless stated otherwise, use $ \mu=1 $, $ m=1 $, $ g=9.81 $, $ l=1 $, $ \theta_0=85^o $ and $ \dot\theta_0=0 $.

Problem 2: Euler's method

Write a Python program that solves the ODE in \eqref{eq:1} with the specified initial conditions using Euler's method (remember to use radians). Experiment with different time steps $ \Delta t $, and friction coefficients $ \mu $. How do they affect the stability of the solution? Carry out the computation for a few cycles, and plot the amplitude as a functions of $ t $.

Solution.

import numpy as np
import matplotlib.pylab as plt
from myanimation import animatePendulum

#### parameters ####

m = 1.
g = 9.81
L = 1.
mu = 1.

#### initial conditions ####
theta_0 = np.pi*85./180
dtheta_0 = 0.

#### arrays for time and solution vectors ####
t_start = 0
t_end = 10.
N = 1000

T = np.linspace(t_start, t_end, N + 1)
dt = T[1] - T[0]
Theta0 = np.zeros_like(T) # array to contain solutions for theta
Theta1 = np.zeros_like(T) # array to contain solutions for d(theta)/dt

Theta0[0] = theta_0 # set initial condition for theta
Theta1[0] = dtheta_0 # set initial condition for d(theta)/dt

#### solving loop for euler's method
for n, t in enumerate(T[:-1]):
    
    f0 = Theta1[n]
    f1 = -mu*Theta1[n]/m - g*np.sin(Theta0[n])/L
    Theta0[n + 1] = Theta0[n] + dt*f0
    Theta1[n + 1] = Theta1[n] + dt*f1

plt.figure()
plt.plot(T, Theta0)
plt.show()

E = 0.5*(Theta1**2 - dtheta_0**2) + np.cos(theta_0) - np.cos(Theta0)

animatePendulum(Theta0, T, l=1)

Problem 3: Heun's method

a) Solve Problem 2 using Heun's method.

b) Write a function that detects the amplitude peaks, and use this to calculate the period of the pendulum. Is it constant?

Solution.

import numpy as np
import matplotlib.pylab as plt
#from myanimation import animatePendulum

def findPeaks(X, V):
    peakList = []
    
    for n, v in enumerate(V[:-1]):
        if np.sign(v) != np.sign(V[ n + 1]):
            peakList.append(X[n])
    
    return np.array(peakList)

#### parameters ####

m = 1.
g = 9.81
L = 2.
mu = 1

#### initial conditions ####
theta_0 = np.pi*85./180
dtheta_0 = 0.

#### arrays for time and solution vectors ####
t_start = 0
t_end = 10.
N = 10000

T = np.linspace(t_start, t_end, N + 1)
dt = T[1] - T[0]
Theta0 = np.zeros_like(T) # array to contain solutions for theta
Theta1 = np.zeros_like(T) # array to contain solutions for d(theta)/dt

Theta0[0] = theta_0 # set initial condition for theta
Theta1[0] = dtheta_0 # set initial condition for d(theta)/dt

#### solving loop for heun's method
for n, t in enumerate(T[:-1]):
    
    f0 = Theta1[n]
    f1 = -mu*Theta1[n]/m - g*np.sin(Theta0[n])/L

    theta0_p = Theta0[n] + dt*f0
    theta1_p = Theta1[n] + dt*f1

    f0_p = theta1_p
    f1_p = -mu*theta1_p/m - g*np.sin(theta0_p)/L

    Theta0[n + 1] = Theta0[n] + dt*(f0 + f0_p)/2.
    Theta1[n + 1] = Theta1[n] + dt*(f1 + f1_p)/2.
    

plt.figure()
plt.plot(T, Theta0)
plt.plot(T, Theta1)
plt.plot(findPeaks(T, Theta1), np.zeros_like(findPeaks(T, Theta1)), 'o')
print findPeaks(T, Theta1)[1:] - findPeaks(T, Theta1)[:-1]
plt.show()

#animatePendulum(Theta0, T, l=1)

Problem 4: (Optional) General second order scheme

a) Solve Problem 2 using Midpoint method and Ralston's method. See: "See http://folk.ntnu.no/fredf/teaching/tkt4140/exercises/pythonExercises/pexercise2/TKT4140%20Lec%204%20RKmethods%20(1).pdf

c) create an animation similar to the one above.

Solution.

import numpy as np
import matplotlib.pylab as plt
#from myanimation import animatePendulum

#### parameters ####

m = 1.
g = 9.81
L = 2.
mu = 1

#### initial conditions ####
theta_0 = np.pi*85./180
dtheta_0 = 0.

#### arrays for time and solution vectors ####
t_start = 0
t_end = 10.
N = 10000

T = np.linspace(t_start, t_end, N + 1)
dt = T[1] - T[0]
Theta0 = np.zeros_like(T) # array to contain solutions for theta
Theta1 = np.zeros_like(T) # array to contain solutions for d(theta)/dt

Theta0[0] = theta_0 # set initial condition for theta
Theta1[0] = dtheta_0 # set initial condition for d(theta)/dt

solvingScheme = 'midpoint'

if solvingScheme == 'midpoint':
    #### solving loop for midpoint's method
    for n, t in enumerate(T[:-1]):
        
        K1_0 = Theta1[n]
        K1_1 = -mu*Theta1[n]/m - g*np.sin(Theta0[n])/L
    
        theta0_p = Theta0[n] + 0.5*dt*K1_0
        theta1_p = Theta1[n] + 0.5*dt*K1_1
    
        K2_0 = theta1_p
        K2_1 = -mu*theta1_p/m - g*np.sin(theta0_p)/L
    
        Theta0[n + 1] = Theta0[n] + dt*(K2_0)
        Theta1[n + 1] = Theta1[n] + dt*(K2_1)

else:
    #### solving loop for Ralston's method
    for n, t in enumerate(T[:-1]):
        
        K1_0 = Theta1[n]
        K1_1 = -mu*Theta1[n]/m - g*np.sin(Theta0[n])/L
    
        theta0_p = Theta0[n] + 0.75*dt*K1_0
        theta1_p = Theta1[n] + 0.75*dt*K1_1
    
        K2_0 = theta1_p
        K2_1 = -mu*theta1_p/m - g*np.sin(theta0_p)/L
    
        Theta0[n + 1] = Theta0[n] + dt*(K1_0/3. +2*K2_0/3.)
        Theta1[n + 1] = Theta1[n] + dt*(K1_1/3. + 2*K2_1/3)
    

plt.figure()
plt.plot(T, Theta0)
plt.plot(T, Theta1)
plt.show()

#animatePendulum(Theta0, T, l=1)

# python_exercises/pexercise2/myanimation.py

import matplotlib.pyplot as plt
from math import sin, cos, pi
from matplotlib import animation

def animatePendulum(theta, t, l=1):
            
    # set up figure and animation
    
    fig = plt.figure()
    ax = fig.add_subplot(111, aspect='equal', autoscale_on=False,
                         xlim=(-1.1, 1.1), ylim=(-1.1, 1.1))
    #ax.grid()
    
    line, = ax.plot([], [], '-o', lw=2)
    time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes)
    
    theta_text =  ax.text(0.02, 0.85, '', transform=ax.transAxes)
    def init():
        """initialize animation"""
        line.set_data([], [])
        time_text.set_text('')
        
        theta_text.set_text('')
        return line, time_text
    
    def animate(i):
        """perform animation step"""
        x = l*sin(theta[i])
        y = -l*cos(theta[i])
        deg = theta[i]*180/pi
        
        line.set_data([0, x], [0, y])
        time_text.set_text('time = %.1f [-]' % t[i])
        theta_text.set_text(r'$\theta$ = %.0f $^o$' % deg)
        return line, time_text
    
    # choose the interval based on dt and the time to animate one step
    Npictures = len(t)
    ani = animation.FuncAnimation(fig, animate,init_func=init,frames=Npictures,interval=1,blit=False)
    

    #Writer = animation.writers['ffmpeg']
    #writer = Writer(fps=30, metadata=dict(artist='Me'), bitrate=1800)
    #ani.save('pendulum.mov', writer=writer)
    
    plt.show()