Multiple Boundary Layers

# Our problem. show(r"$\epsilon^2 y'' + \epsilon x y' - y = - e^x, \quad y(0) = 2, y(1) = 1.$") 

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{$\epsilon^2 y'' + \epsilon x y' - y = - e^x, \quad y(0) = 2, y(1) = 1.$}

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{$\epsilon^2 y'' + \epsilon x y' - y = - e^x, \quad y(0) = 2, y(1) = 1.$}

# Perturbation solution. (x, epsilon) = var('x, epsilon') y_perturbation0 = exp(x) + exp(- x / epsilon) - (exp(1) - 1) * exp(- (sqrt(5) - 1) / 2 * (1 - x) / epsilon) show(r'$y_0(x) = %s.$' % latex(y_perturbation0)) 

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{$y_0(x) = -{\left(e - 1\right)} e^{\left(\frac{{\left(\sqrt{5} - 1\right)} {\left(x - 1\right)}}{2 \, \epsilon}\right)} + e^{x} + e^{\left(-\frac{x}{\epsilon}\right)}.$}

\newcommand{\Bold}[1]{\mathbf{#1}}\hbox{$y_0(x) = -{\left(e - 1\right)} e^{\left(\frac{{\left(\sqrt{5} - 1\right)} {\left(x - 1\right)}}{2 \, \epsilon}\right)} + e^{x} + e^{\left(-\frac{x}{\epsilon}\right)}.$}

# Numerical solution # # Writing the 2nd-order equation as the 1st-order system: # y' = yp # yp' = - x * yp / epsilon + y / epsilon^2 - exp(x) / epsilon^2. def RHS(Y, x, epsilon): 'Right-hand side of the ODE system.' (y, yp) = Y dy = yp dyp = (- x * yp + (y - exp(x)) / epsilon) / epsilon return (dy, dyp) # Boundary conditions. xl = 0. yl = 2. xr = 1. yr = 1. 

# Boundary-value problem solver. import scipy.integrate import scipy.optimize class BVP: ''' Solve the 2nd-order boundary-value problem y'' = f(y, y', x), y(xl) = yl, y(xr) = yr, using a shooting method. RHS(Y, x, *args) is a standard ODE function for scipy.integrate.odeint. x is the points at which to solve, (xl, ..., xr). yl is the left-hand boundary condition. yr is the right-hand boundary condition. Since ODE solvers solve initial-value problems y'' = f(y, y', x), y(xl) = yl, y'(xl) = ypl, we will find the value for ypl that gives the solution that satisfies the right-hand boundary condition y(xr) = yr. ''' def __init__(self, RHS, x, yl, yr, *args): self.RHS = RHS self.x = x self.yl = yl self.yr = yr self.args = args def _solve_IVP(self, ypl_guess): Y = scipy.integrate.odeint(self.RHS, (self.yl, ypl_guess), self.x, args = self.args) (y, yp) = zip(*Y) return y def _error_in_yr(self, ypl_guess): y = self._solve_IVP(ypl_guess) return y[-1] - self.yr def solve(self, ypl_initialguess = 0.): # Use a Newton-type method. ypl = scipy.optimize.newton(self._error_in_yr, ypl_initialguess) # Return the solution. return self._solve_IVP(ypl) 

# Interactively plot solution & approximation. # Number of points on each curve. npoints = 1001 # Plotting styles. perturbation_style = {'color': 'red', 'linestyle': '-', 'alpha': 0.5} numerical_style = {'color': 'black', 'linestyle': '--', 'alpha': 0.5} endpoints_style = {'color': 'black', 'size': 20} # Make the plots. @interact def plot_solutions(log10epsilon = slider(log(0.03, 10), log(0.3, 10), default = -1, label = '$\log_{10} \epsilon$')): epsilon = 10^log10epsilon html(r'$$\epsilon = %s$$' % latex(epsilon)) plots = [] # Perturbation solution. plots.append(plot(y_perturbation0(epsilon = epsilon), x, 0, 1, plot_points = npoints, **perturbation_style)) # Numerical solution. x_numerical = srange(xl, xr, (xr - xl) / (npoints - 1), include_endpoint = true) y_numerical = BVP(RHS, x_numerical, yl, yr, epsilon).solve() plots.append(line(zip(x_numerical, y_numerical), **numerical_style)) # Add boundary points. plots.append(point([(xl, yl), (xr, yr)], **endpoints_style)) show(sum(plots), axes_labels = (r'$x$', r'$y(x)$'))