BVP layer -- Solve a singular perturbation problem¶
AMath 585, Winter Quarter 2020 at the University of Washington. Developed by R.J. LeVeque and distributed under the BSD license. You are free to modify and use as you please, with attribution.
These notebooks are all available on Github.
Solve the nonlinear BVP $$ \epsilon u''(x) + u(x)(u'(x) - 1) = f(x) $$ with Dirichlet boundary conditions.
In this notebook we consider the case where $\epsilon > 0$ is very small, the singular perturbation problem discussed in Section 2.17 of the textbook.
Continuation in both $\epsilon$ and the size of the system are illustrated as well.
%matplotlib inline
from pylab import *
from scipy.interpolate import interp1d # used for continuation
The BVP_nonlinear.py module contains the function solve_bvp_nonlinear that is illustrated in the notebook BVP_nonlinear.html.
from BVP_nonlinear import solve_bvp_nonlinear
ax = 0.; alpha = -1.; ainfo = (ax, alpha)
bx = 1.; beta = 1.5; binfo = (bx, beta)
f = lambda x: zeros(x.shape) # the zero function
xbar = 0.5*(ax+bx-beta-alpha)
print('For small epsilon we expect a layer near xbar = %.2f' % xbar)
epsilon = 0.1
m = 49
u0_func = lambda x: alpha + x * (beta-alpha) / (bx-ax)
x,u = solve_bvp_nonlinear(epsilon, f, ainfo, binfo, m, u0_func)
For smaller epsilon we see slower convergence and some wonky looking iterates, but the final solution looks smooth with a layer where we expect it:
epsilon = 0.01
m = 199
u0_func = lambda x: alpha + x * (beta-alpha) / (bx-ax)
x,u = solve_bvp_nonlinear(epsilon, f, ainfo, binfo, m, u0_func)
figure()
plot(x,u,'rx-')
title('Final approximate solution')
Note that we have many grid points in the layer. If we tried this value of epsilon with fewer grid points the final solution won't look so nice, and/or Newton's method might not converge.
Here we see that it doesn't converge when the starting guess is the linear function:
epsilon = 0.01
m = 49
u0_func = lambda x: alpha + x * (beta-alpha) / (bx-ax)
x,u = solve_bvp_nonlinear(epsilon, f, ainfo, binfo, m,
u0_func, max_iter=15)
epsilon = 0.05
m = 49
u0_func = lambda x: alpha + x * (beta-alpha) / (bx-ax)
x,u = solve_bvp_nonlinear(epsilon, f, ainfo, binfo, m, u0_func)
u_05 = u
Now we can reduce epsilon and use the result just computed as the initial guess instead of a linear function.
Note that u0_func is a function that just returns this vector regardless of what x is passed in, so this wouldn't work if we wanted to increase m (which we will do later).
epsilon = 0.01
m = 49
u0_func = lambda x: u_05
x,u = solve_bvp_nonlinear(epsilon, f, ainfo, binfo, m, u0_func)
u_01 = u
figure()
plot(x,u,'rx-')
title('Final approximate solution')
Now Newton's method converges and we get something reasonable, even though there are only a few points in the interior layer.
We saved the solution above, which we can now use as initial guess for an even smaller epsilon:
epsilon = 0.007
m = 49
u0_func = lambda x: u_01
x,u = solve_bvp_nonlinear(epsilon, f, ainfo, binfo, m, u0_func, max_iter=20)
figure()
plot(x,u,'rx-')
title('Final approximate solution')
Continuation in m (refining the grid)¶
Newton's method wasn't converging in the test just done, and clearly we don't have enough grid points in the layer for this value of epsilon. So we might want to increase m. To do so we can define a function based on one of our previous converged solutions that we can then pass in as u0_func and that can be evaluated on a new finer grid.
To do so we use the scipy.interpolate.interp1d function that was imported at the top of this notebook, to define a piecewise linear function that interpolates the values (x[i], u_01[i]):
u_01_func = interp1d(x,u_01,'linear')
epsilon = 0.007
m = 99
x,u = solve_bvp_nonlinear(epsilon, f, ainfo, binfo, m, u0_func=u_01_func)
figure()
plot(x,u,'rx-')
title('Final approximate solution')
With this approach we can go to smaller epsilon, on a suitable fine grid:
epsilon = 0.001
m = 499
x,u = solve_bvp_nonlinear(epsilon, f, ainfo, binfo, m, u0_func=u_01_func)
figure()
plot(x,u,'rx-')
title('Final approximate solution')
