Ir al contenido principal

# Numerical methods challenge: Day 20

During October (2017) I will write a program per day for some well-known numerical methods in both Python and Julia. It is intended to be an exercise then don't expect the code to be good enough for real use. Also, I should mention that I have almost no experience with Julia, so it probably won't be idiomatic Julia but more Python-like Julia.

## Shooting method

Today we have the shooting method. This is a method for solving boundary value problems by turning them into a sequence of initial value problems.

Loosely speaking, for a second order equation we have

\begin{equation*} y''(x) = f(x, y(x), y'(x)),\quad y(x_0) = y_0,\quad y(x_1) = y_1\, , \end{equation*}

and we solve the sequence of problems

\begin{equation*} y''(x) = f(x, y(x), y'(x)),\quad y(x_0) = y_0,\quad y'(x_0) = a \end{equation*}

until the function $F(a) = y(x_1; a) - y_1$ has a root.

We will test the method with the boundary value problem

\begin{equation*} y''(x) = \frac{3}{2} y^2,\quad w(0) = 4,\quad w(1) = 1\, . \end{equation*}

Following are the codes.

### Python

from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import newton
from scipy.integrate import odeint

def shooting(dydx, x, x0, xf, shoot=None):
if shoot is None:
shoot = np.random.uniform(-20, 20)
F = lambda s, x0, xf, x: odeint(dydx, [x0, s], x)[-1, 0] - xf
shoot = newton(F, shoot, args=(x0, xf, x))
y = odeint(dydx, [x0, shoot], x)
return y[:, 0], shoot

func = lambda y, t: [y, 1.5*y**2]
x = np.linspace(0, 1, 1000)
y, shoot = shooting(func, x, 4, 1, shoot=-5)
plt.plot(x, y)
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.show()


### Julia

using PyPlot, DifferentialEquations, Roots

function shooting(dydx, x, y0, yf; shoot=nothing)
if shoot == nothing
shoot = rand(-20:0.1:20)
end
function F(s)
prob = ODEProblem(dydx, [y0, s], (x, x[end]))
sol = solve(prob)
return sol(x[end]) - yf
end
shoot = fzero(F, shoot)
prob = ODEProblem(dydx, [y0, shoot], (x, x[end]))
sol = solve(prob, solveat=x)
return sol(x)[1, :], shoot
end

func(x, y) = [y, 1.5*y^2]
x = linspace(0, 1, 1000)
y, shoot = shooting(func, x , 4.0, 1.0, shoot=-5.0)
plot(x, y)


In both cases the result is the following plot, and the slope is -7.9999999657800833. We should mention that the convergence of the method relies on the selection of initial guesses. For example, if we choose as initial parameter -50 in the previous problem, we obtain a completely differente solution.

y, shoot = shooting(func, x , 4.0, 1.0, shoot=-50.0) And the obtained slope is -35.858547970130971.

### Comparison Python/Julia

Regarding number of lines we have: 20 in Python and 23 in Julia. The comparison in execution time is done with %timeit magic command in IPython and @benchmark in Julia.

For Python:

%timeit shooting(func, x, 4, 1, shoot=-5)


with result

100 loops, best of 3: 1.9 ms per loop


For Julia:

@benchmark shooting(func, x, 4.0, 1.0, shoot=-5.0)


with result

BenchmarkTools.Trial:
memory estimate:  4.18 MiB
allocs estimate:  78552
--------------
minimum time:     10.065 ms (0.00% GC)
median time:      10.593 ms (0.00% GC)
mean time:        11.769 ms (9.28% GC)
maximum time:     22.252 ms (48.58% GC)
--------------
samples:          425
evals/sample:     1


In this case, we can say that the Python code is roughly 5 times faster than Julia. Although, the codes are more different than in the other challenges. For example, I am not using newton in Julia.