jeffery_hamel.py

"""
This script solves the 3rd-order Jeffery-Hamel ODE:

f''' + 2*Re*alpha*f*f' + 4*alpha*alpha*f' = 0
f(0) = 1
f(1) = 0
f'(0) = 0

as three ODEs:

y0' = y1
y1' = y2
y2' = -2 * Re * alpha * y0 * y1 - 4 * alpha**2 * y1

y0(0) = 1
y0(1) = 0
y1(0) = 0

Installation:
* Use pip to install the scikits.bvp_solver package,
  sudo -E pip install scikits.bvp_solver

Running:
* python jeffery_hamel.py
"""
import scikits.bvp_solver
import numpy as np


# Set global constants which are used below.
Re = -80
alpha = 5 * np.pi/180


def function(X, Y):
    """
    'function' argument to bvp_solver.ProblemDefinition().
    """
    return np.array([Y[1],                                                 # y0'
                     Y[2],                                                 # y1'
                     -2 * Re * alpha * Y[0] * Y[1] - 4 * alpha**2 * Y[1]]) # y2'

def function_derivative(X, Y):
    """
    'function_derivative' argument to bvp_solver.ProblemDefinition().
    """
    dFdY = np.array([[0.0,                    1.0,                                   0.0],
                     [0.0,                    0.0,                                   1.0],
                     [-2 * Re * alpha * Y[1], -2 * Re * alpha * Y[0] - 4 * alpha**2, 0.0]])
    return dFdY

def boundary_conditions (Ya, Yb):
    """
    'boundary_conditions' argument to bvp_solver.ProblemDefinition().
    """
    BCa = np.array([Ya[0] - 1.0,    # y0(a) = 1
                    Ya[1] - 0.0])   # y1(a) = 0
    BCb = np.array([Yb[0] - 0.0])   # y0(b) = 0
    return BCa, BCb


def boundary_conditions_derivative(YA, YB):
    """
    'boundary_conditions_derivative' argument to bvp_solver.ProblemDefinition().
    """
    dBCdYa = np.array([[1.0, 0.0, 0.0],   # left BC 0
                       [0.0, 1.0, 0.0]])  # left BC 1
    dBCdYb = np.array([[1.0, 0.0, 0.0]])  # right BC 0
    return dBCdYa, dBCdYb


def guess(X):
    """
    'solution_guess' argument to bvp_solver.solve().
    """
    return np.array([1-X, -1., 0.])


# Set up the problem
problem = scikits.bvp_solver.ProblemDefinition(num_ODE = 3,
                                               num_parameters = 0,
                                               num_left_boundary_conditions = 2,
                                               boundary_points = (0, 1),
                                               function = function,
                                               boundary_conditions = boundary_conditions,
                                               function_derivative = function_derivative,
                                               boundary_conditions_derivative = boundary_conditions_derivative)

# Solve the problem.
# We pass 'trace=2' (argument can be in the range [0,2]) to solve() to
# control the amount of output from the algorithm during the solve.
solution = scikits.bvp_solver.solve(problem,
                                    solution_guess = guess,
                                    tolerance = 1.0e-6,
                                    max_subintervals = 5000,
                                    trace = 2)

# The solution object can reconstruct the solution at any desired points.
# y is a numpy ndarray.  "11" points corresponds to 10 elements.
x = np.linspace(problem.boundary_points[0], problem.boundary_points[1], 11)
y = solution(x)

# Print the solution and first derivative values at the requested points.
for j in xrange(y.shape[1]):
    print('{:.16e}, {:.16e},'.format(y[0,j], y[1,j]))