Python Exercise 3

Fredrik Eikeland Fossan

Mar 11, 2019

Shooting technique and discretization methods for solving ODE's

Problem 1: Shooting technique

Given the following boundary value problem $$ \begin{equation} \label{eq:nonlin1} \begin{array}{l} y''= -2 \,x\, \left(y'\right)^2 \\ y(0) = 2 \, , \quad \ y(0.9) = \frac{1}{2}\left[ln \, \left( 0.1\right) - ln \, \left( 1.9 \right) \right] +2 \approx 0.52778 \end{array} \end{equation} $$

where $ y=y(x) $. The analytical solution is: $$ \begin{equation} \label{eq:sol1} y = \frac{1}{2}\left(ln \, \lvert x - 1\rvert - ln \, \lvert x + 1\rvert \right) +2 \end{equation} $$

a) Solve the boundary value problem given in Eq. \eqref{eq:nonlin1}. Choose $ S_0 = -0.5 $, and $ S_1 = -0.6 $, where $ S = y'(0) $. Use the secant method to find $ S_2 $, $ S_3 $ .... Choose for instance rk4 solver, and segmenting $ x $ from $ 0 $ to $ 0.9 $ with 91 points.

Hint 1.

Importing matplotlib in the following manner, and adding the line below will make your figures pop up "in front" of the Liclipse window:

import matplotlib; matplotlib.use('Qt4Agg')
import matplotlib.pylab as plt
plt.get_current_fig_manager().window.raise_()

The analytical solution can be created by:

N = 90
x = np.linspace(0, 0.9, N + 1)

y_analytic = 0.5*(np.log(np.abs(x-1))-np.log(np.abs(x+1))) + 2

Hint 2.

After solving with 7 different $ S $ values, $ S_6 \approx -0.998 $, $ \phi_6 \approx 0.005 $. $ S^* = -1 $, is the correct value, $ \phi\left(S^*\right)=0 $.

Figure 1: Curves from the different iterations.

b) Plot the $ \phi $ function from $ S=-0.5 $ to $ S=-1.23 $. Segment with 51 points.

Hint.

The following template script could be used to plot the $ \phi $ function

s_start, s_end = -0.5, -1.23
sList = np.linspace(s_start, s_end, 51)
phiList = np.zeros_like(sList)

for n, s in enumerate(sList):
    Y_0 = 
    Y_shoot = 
    y_shoot = # extract y, and not y'
    phiList[n] = 
    
plt.figure()
plt.plot(sList, phiList)
plt.plot(sList, np.zeros_like(sList), 'k--')
plt.xlim([s_start, s_end])
plt.ylim([np.min(phiList), np.max(phiList)])
plt.xlabel('s')
plt.ylabel(r'$\phi$')

Figure 2: Plot off $ \phi $ function.

Solution.

# python_exercises/pexercise4/problem1.py;ODEschemes.py @ git@lrhgit/tkt4140/src/src-ch2/ODEschemes.py;
#import matplotlib; matplotlib.use('Qt4Agg')
import matplotlib.pylab as plt
#plt.get_current_fig_manager().window.raise_()

import numpy as np
import math
from ODEschemes import rk4

LNWDT=3; FNT=15
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT


def dUfunc(Y, x):
    y0 = Y[0]
    y1 = Y[1]
    return np.array([y1, -2*x*y1**2])

N = 90
x = np.linspace(0, 0.9, N + 1)

y_analytic = 0.5*(np.log(np.abs(x-1))-np.log(np.abs(x+1))) + 2

beta = 0.52778 #0.5*(math.log(0.1) - math.log(1.9)) + 2

s0 = -0.5
Y_0 = [2, s0]

Y_shoot = rk4(dUfunc, Y_0, x)
y_shoot = Y_shoot[:, 0]
phi0 = y_shoot[-1] - beta

s1 = -0.6
plt.figure()
grey = '0.75'
linestyles = [grey, 'r--', 'g--', 'b--', 'y--', 'm--', 'c--']

plt.plot(x, y_analytic, 'k')
plt.plot(x, y_shoot, linestyles[0])
legendList = ['analytic', 'it = 0']
for n in range(6):
    Y_0 = [2., s1]
    Y_shoot = rk4(dUfunc, Y_0, x)
    print Y_shoot[-1, 1]
    y_shoot = Y_shoot[:, 0]
    phi1 = y_shoot[-1] - beta
    print "s1: {0}; phi1: {1}".format(s1, phi1)
    ds = -phi1 *(s1 - s0)/float(phi1 - phi0)
    s0 = s1
    s1 = s1 + ds
    phi0 = phi1
    plt.plot(x, y_shoot, linestyles[n + 1])
    legendList.append('it = {0}'.format(n + 1))

plt.xlabel('x')
plt.ylabel('y')
plt.legend(legendList, frameon=False, loc='best')
#plt.savefig('fig/Problem1_it.png', transparent=True)

plt.figure()
plt.plot(x, y_analytic, 'k')
plt.plot(x, y_shoot, 'r--')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(['analytic', 'numeric'], frameon=False, loc='best')

#plt.savefig('fig/Problem1_final.png', transparent=True)

#### plot Phi-funcion ####
s_start, s_end = -0.5, -1.23
sList = np.linspace(s_start, s_end, 51)
phiList = np.zeros_like(sList)

for n, s in enumerate(sList):
    Y_0 = [2, s]
    Y_shoot = rk4(dUfunc, Y_0, x)
    y_shoot = Y_shoot[:, 0] # extract y, and not y'
    phiList[n] = y_shoot[-1] - beta
    
plt.figure()
plt.plot(sList, phiList)
plt.plot(sList, np.zeros_like(sList), 'k--')
plt.xlim([s_start, s_end])
plt.ylim([np.min(phiList), np.max(phiList)])
plt.xlabel('s')
plt.ylabel(r'$\phi$')
#plt.savefig('fig/Problem1_phi.png', transparent=True)
plt.show()

Problem 2: Solving ODE's using Discretization methods

Solve problem 1 c) from the last theory exercise. Solve with $ h=0.25 $ (compare with your answers from the theory exercise), $ h=0.1 $ and $ h=0.01 $. Add the boundaries and plot the solutions.

Hint 1.

The function scipy.sparse.diags is very convenient to use when creating sparse matrices. Consider the following matrix: $$ \begin{equation} A = \left[\begin{array}{ccccc} 2& 3 & 0 & 0 & 0 \\ 1& 2 & 3 & 0 & 0 \\ 0& 1 & 2 & 3 & 0 \\ 0& 0 & 1 & 2 & 3\\ 0& 0 & 0 & 1 & 2 \end{array}\right] \label{_auto1} \end{equation} $$

This could be created by

import numpy as np
import scipy
import scipy.sparse
import scipy.sparse.linalg
N = 5
A = scipy.sparse.diags([1, 2, 3], [-1, 0, 1], shape=(N, N), format='csc') # format='csc' works well with scipy.sparse.linalg.spsolve
A_numpy = A.toarray() # convert to numpy array

or alternatively:

N = 5
subDiag = np.ones(N - 1)
mainDiag = 2*np.ones(N)
superDiag = 3*np.ones(N - 1)
A = scipy.sparse.diags([subDiag, mainDiag, superDiag], [-1, 0, 1], format='csc')

To solve an equation $ A \cdot x = d $ where $ x = [x1, x2, x3, x4, x5]^T $ and $ d =[ 8., 14., 20., 26., 14.]^T $:

d = np.array([  8,  14,  20,  26,  14])

x = scipy.sparse.linalg.spsolve(A, d) # sparse solver.
x2 = np.linalg.solve(A_numpy, d) # numpy linalg solver

which gives $ x = [1, 2, 3, 4, 5]^T $

Hint 2.

Figure 3: Solution with $ h=0.1 $.

Solution.

# python_exercises/pexercise4/problem2.py;

#import matplotlib; matplotlib.use('Qt4Agg')
import matplotlib.pylab as plt
#plt.get_current_fig_manager().window.raise_()
import numpy as np
import scipy
import scipy.sparse
import scipy.sparse.linalg

LNWDT=3; FNT=15
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT


N = 3
h = 1./(N + 1)
Beta2 = 4.

x = np.linspace(0, 1, N + 2)


i_list = np.linspace(0, N + 1, N + 2)


MainDiag = -(2*i_list[1:-1] + Beta2*h)
subDiag = i_list[2:-1] - 1
supDiag = i_list[1:-2] + 1
d = np.zeros(N)
d[-1] = -(N + 1)

A = scipy.sparse.diags([subDiag, MainDiag, supDiag], [-1, 0, 1], format='csc')

ThetaSol = scipy.sparse.linalg.spsolve(A, d)

theta1 = ThetaSol[0]
theta0 = theta1/(1 + (Beta2/2)*h + (Beta2**2/12)*h**2)
thetaEnd = 1
Theta = np.append(theta0, ThetaSol)
Theta = np.append(Theta, thetaEnd)

print Theta
plt.figure()
plt.plot(x, Theta)
plt.xlabel('x')
plt.ylabel(r'$\theta$')
#plt.savefig('fig/Problem2.png', transparent=True)
plt.show()

Problem 3: Linearization

solve Problem 1 by discretizing both $ y'' $ and $ y' $ using central differences.

a) Linearize the term $ \left( y^{m + 1}_{i + 1} - y^{m + 1}_{i - 1}\right)^2 $ by $ \left( y^{m + 1}_{i + 1} - y^{m + 1}_{i - 1}\right)^2 \approx \left( y^{m}_{i + 1} - y^{m}_{i - 1}\right)^2 $, where $ m $ denotes the iteration. Note that this is slightly different than "metoden med etterslep", and the matrix should now only consist of constants, whereas the RHS will contain $ y $-values from the previous iteration. Choose for instance $ h=0.1 $, which should give 8 equations and 8 unknowns.

Hint 1.

To start the iteration we need to choose $ Y^0 $ (m = 0). A straight line between the boundary values could be a reasonable starting point:

N = 8 # number of unknowns
x = np.linspace(0, 0.9, N + 2)
h = 0.9/(N + 1)
y_analytic = 0.5*(np.log(np.abs(x-1))-np.log(np.abs(x+1))) + 2

y_0, y_End = y_analytic[0], y_analytic[-1] # boundaries

Y0 = np.linspace(y_0, y_End, N + 2) # first value of Y

Hint 2.

Figure 4: Solution of a) after two iterations, with $ h=0.1 $.

Solution.

# python_exercises/pexercise4/problem3a.py;

#import matplotlib; matplotlib.use('Qt4Agg')
import matplotlib.pylab as plt
#plt.get_current_fig_manager().window.raise_()
import numpy as np
import scipy.linalg
import scipy.sparse
import scipy.sparse.linalg

LNWDT=2; FNT=15
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT


N = 8
x = np.linspace(0, 0.9, N + 2)
h = 0.9/(N + 1)

y_analytic = 0.5*(np.log(np.abs(x-1))-np.log(np.abs(x+1))) + 2

x_unknown = x[1:-1] # all x-values except BC's

y0, yEnd = y_analytic[0], y_analytic[-1]

Y0 = np.linspace(y0, yEnd, N + 2) # initial guess

for n in range(2):
    Y0Plus = Y0[2:]
    Y0Minus = Y0[0:-2]
    d = -0.5*x_unknown*(Y0Plus-Y0Minus)**2
    d[0]= d[0] - y0
    d[-1]= d[-1] - yEnd
    
    A = scipy.sparse.diags([1, -2, 1], [-1, 0, 1], shape=(N, N), format='csc')
    
    Y1 = scipy.sparse.linalg.spsolve(A, d)
    
    Y1Full = np.append(y0, Y1)
    Y1Full = np.append(Y1Full, yEnd)
    
    Y0 = Y1Full
    

plt.figure()
plt.plot(x, y_analytic, 'k')
plt.plot(x, Y1Full, 'r--')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(['analytic', 'numeric'], frameon=False)
#plt.savefig('fig/Problem3a.png', transparent=True)
plt.show()

b) Try also to Linearize the term $ \left( y^{m + 1}_{i + 1} - y^{m + 1}_{i - 1}\right)^2 $ by using Newton linearization. Solve the linearized system, and compare with a). Newton linearization should converge faster.

Hint 1.

After employing Newton linearization and performing some algebraic manipulation: $$ \begin{equation} \label{eq:newtonLin} \left( y^{m + 1}_{i + 1} - y^{m + 1}_{i - 1}\right)^2 \approx 2 \alpha_i^m \, y^{m + 1}_{i + 1} - 2 \alpha_i^m \, y^{m + 1}_{i - 1} - \left(\alpha_i^m \right)^2 \end{equation} $$

where $ \alpha_i^m = y^{m}_{i + 1} - y^{m}_{i - 1} $. Inserting to the discretized equation should give: $$ \begin{equation} \label{eq:EQnewtonLin} \left(1-x_i \,\alpha_i^m\right)y^{m + 1}_{i - 1} - 2y^{m + 1}_{i} + \left(1+x_i \,\alpha_i^m\right)y^{m + 1}_{i + 1}=\frac{1}{2}x_i \,\left(\alpha_i^m\right)^2 \end{equation} $$ Try to show this yourself.

Hint 2.

Figure 5: Solution of a) after two iterations, with $ h=0.1 $.

Solution.

#import matplotlib; matplotlib.use('Qt4Agg')
import matplotlib.pylab as plt
#plt.get_current_fig_manager().window.raise_()
import numpy as np
import scipy.linalg
import scipy.sparse
import scipy.sparse.linalg

LNWDT=2; FNT=15
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT

N = 8
x = np.linspace(0, 0.9, N + 2)
h = 0.9/(N + 1)

y_analytic = 0.5*(np.log(np.abs(x-1))-np.log(np.abs(x+1))) + 2

x_unknown = x[1:-1]

y0, yEnd = y_analytic[0], y_analytic[-1]

Y0 = np.linspace(y0, yEnd, N + 2)

for n in range(2):
    Y0Plus = Y0[2:]
    Y0Minus = Y0[0:-2]
    alpha = Y0Plus - Y0Minus
    
    MainDiag = -2*np.ones(N)
    supDiag = 1 + x_unknown[:-1]*alpha[:-1]
    subDiag = 1 - x_unknown[1:]*alpha[1:]
    d = 0.5*x_unknown*(alpha)**2
    d[0]= d[0] - y0*(1 - x_unknown[0]*alpha[0])
    d[-1]= d[-1] - yEnd*(1 + x_unknown[-1]*alpha[-1])
    
    A = scipy.sparse.diags([subDiag, MainDiag, supDiag], [-1, 0, 1], format='csc')
    
    Y1 = scipy.sparse.linalg.spsolve(A, d)
    
    Y1Full = np.append(y0, Y1)
    Y1Full = np.append(Y1Full, yEnd)
    
    Y0 = Y1Full
    
print len(x), len(Y0), len(Y1Full)
plt.figure()
plt.plot(x, y_analytic, 'k')
plt.plot(x, Y1Full, 'r--')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(['analytic', 'numeric'], frameon=False)
#plt.savefig('fig/Problem3b.png', transparent=True)
plt.show()