Python Exercise 3
Feb 26, 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.
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
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.
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.
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.
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 \)
Figure 3: Solution with \( h=0.1 \).
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.
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
Figure 4: Solution of a) after two iterations, with \( h=0.1 \).
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.
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.
Figure 5: Solution of a) after two iterations, with \( h=0.1 \).
