%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sn
import pycollocation
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)
def initial_mesh(yL, yH, num, problem):
ys = np.linspace(yL, yH, num=num)
ns = problem.params['nL'] + np.sqrt(ys)
return ys, ns
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)
solution.result.success
True
n_soln, = solution.evaluate_solution(ys)
plt.plot(ys, n_soln)
plt.show()
n_resids, = solution.evaluate_residual(ys)
plt.plot(ys, n_resids)
plt.show()
n_normalized_resids, = solution.normalize_residuals(ys)
plt.plot(ys, np.abs(n_normalized_resids))
plt.yscale('log')
plt.show()
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
plt.plot(ys, spence_analytic_solution(ys, **params))
plt.plot(ys, n_soln)
plt.show()