Skip to main content

Nicolás' blog

Numerical methods challenge: Day 20

Nicolás Guarín-Zapata

Comments

Also available in:

Español

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], 1.5*y[0]**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[1], x[end]))
        sol = solve(prob)
        return sol(x[end])[1] - yf
    end
    shoot = fzero(F, shoot)
    prob = ODEProblem(dydx, [y0, shoot], (x[1], x[end]))
    sol = solve(prob, solveat=x)
    return sol(x)[1, :], shoot
end


func(x, y) = [y[2], 1.5*y[1]^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.

Solution of the differential equation that satisfy the boundary conditions.

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)
Solution of the differential equation that satisfy the boundary conditions.

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.

Comments

Comments powered by Disqus