``````

In [1]:

%matplotlib inline

``````
``````

In [11]:

import matplotlib.pyplot as plt
import numpy as np
import seaborn as sn

import pycollocation

``````
``````

In [101]:

def spence_model(y, n, alpha, **params):
return [(n**-1 - alpha * n * y**(alpha - 1)) / y**alpha]

def initial_condition(y, n, nL, alpha, **params):
return [n - nL]

def yL(nL, alpha):
return (nL**2 * alpha)**(1 / (1 - alpha))

params = {'nL': 1.0, 'alpha': 0.15}

classic_spence_ivp = pycollocation.problems.IVP(initial_condition, 1, 1, params,
spence_model)

``````
``````

In [102]:

def initial_mesh(yL, yH, num, problem):
ys = np.linspace(yL, yH, num=num)
ns = problem.params['nL'] + np.sqrt(ys)
return ys, ns

``````
``````

In [103]:

bspline_basis = pycollocation.basis_functions.BSplineBasis()
solver = pycollocation.solvers.Solver(bspline_basis)

boundary_points = (yL(**params), 10)
ys, ns = initial_mesh(*boundary_points, num=250, problem=classic_spence_ivp)

tck, u = bspline_basis.fit([ns], u=ys, k=5, s=0)
knots, coefs, k = tck
initial_coefs = np.hstack(coefs)

basis_kwargs = {'knots': knots, 'degree': k, 'ext': 2}
nodes = np.linspace(*boundary_points, num=249)
solution = solver.solve(basis_kwargs, boundary_points, initial_coefs,
nodes, classic_spence_ivp)

``````
``````

In [104]:

solution.result.success

``````
``````

Out[104]:

True

``````
``````

In [106]:

n_soln, = solution.evaluate_solution(ys)
plt.plot(ys, n_soln)
plt.show()

``````
``````

``````
``````

In [107]:

n_resids, = solution.evaluate_residual(ys)
plt.plot(ys, n_resids)
plt.show()

``````
``````

``````
``````

In [108]:

n_normalized_resids, = solution.normalize_residuals(ys)
plt.plot(ys, np.abs(n_normalized_resids))
plt.yscale('log')
plt.show()

``````
``````

``````
``````

In [109]:

def spence_analytic_solution(y, nL, alpha):
"""
Analytic solution to the differential equation describing the signaling
equilbrium of the Spence (1974) model.

"""
# compute analytic solution for worker type given y
D = ((1 + alpha) / 2) * (nL / yL(nL, alpha)**-alpha)**2 - yL(nL, alpha)**(1 + alpha)
return y**(-alpha) * (2 * (y**(1 + alpha) + D) / (1 + alpha))**0.5

``````
``````

In [110]:

plt.plot(ys, spence_analytic_solution(ys, **params))
plt.plot(ys, n_soln)
plt.show()

``````
``````

``````
``````

In [ ]:

``````