Singular Perturbations

(x, xi, epsilon) = var('x, xi, epsilon') 

# Outer solution (Away from x = 0) y0(x) = exp(1 - x) y1(x) = y0(x) + epsilon * (1 - x) / 2 * exp(1 - x) # Display the above html(r'1-Term Outer Solution: $y_0(x) = %s$' % latex(y0(x))) html(r'2-Term Outer Solution: $y_1(x) = %s$' % latex(y1(x))) 

1-Term Outer Solution: y_0(x) = e^{\left(-x + 1\right)}
2-Term Outer Solution: y_1(x) = -\frac{1}{2} \, {\left(x - 1\right)} \epsilon e^{\left(-x + 1\right)} + e^{\left(-x + 1\right)}

1-Term Outer Solution: y_0(x) = e^{\left(-x + 1\right)}
2-Term Outer Solution: y_1(x) = -\frac{1}{2} \, {\left(x - 1\right)} \epsilon e^{\left(-x + 1\right)} + e^{\left(-x + 1\right)}

# Inner solution (Near x = 0; xi = x / epsilon) Y0(xi) = exp(1) - exp(1 - 2 * xi) Y1(xi) = Y0(xi) + epsilon * (1 / 2 * (exp(1) - exp(1 - 2 * xi)) - 2 * xi * exp(1 - 2 * xi)) # Display the above html(r'1-Term Inner Solution: $Y_0(\xi) = %s$' % latex(Y0(xi))) html(r'2-Term Inner Solution: $Y_1(\xi) = %s$' % latex(Y1(xi))) 

1-Term Inner Solution: Y_0(\xi) = -e^{\left(-2 \, \xi + 1\right)} + e
2-Term Inner Solution: Y_1(\xi) = -\frac{1}{2} \, {\left(4 \, \xi e^{\left(-2 \, \xi + 1\right)} + e^{\left(-2 \, \xi + 1\right)} - e\right)} \epsilon - e^{\left(-2 \, \xi + 1\right)} + e

1-Term Inner Solution: Y_0(\xi) = -e^{\left(-2 \, \xi + 1\right)} + e
2-Term Inner Solution: Y_1(\xi) = -\frac{1}{2} \, {\left(4 \, \xi e^{\left(-2 \, \xi + 1\right)} + e^{\left(-2 \, \xi + 1\right)} - e\right)} \epsilon - e^{\left(-2 \, \xi + 1\right)} + e

# Composite solution yC0(x) = y0(x) + Y0(x / epsilon) - y0(0) yC1(x) = y1(x) + Y1(x / epsilon) - y1(0) # Display the above html(r'1-Term Composite Solution: $y_{C 0}(x) = %s$' % latex(yC0(x))) html(r'2-Term Composite Solution: $y_{C 1}(x) = %s$' % latex(yC1(x))) 

1-Term Composite Solution: y_{C 0}(x) = -e^{\left(\frac{-2 \, x}{\epsilon} + 1\right)} + e^{\left(-x + 1\right)}
2-Term Composite Solution: y_{C 1}(x) = -\frac{1}{2} \, {\left(x - 1\right)} \epsilon e^{\left(-x + 1\right)} - \frac{1}{2} \, {\left(\frac{4 \, x e^{\left(\frac{-2 \, x}{\epsilon} + 1\right)}}{\epsilon} + e^{\left(\frac{-2 \, x}{\epsilon} + 1\right)} - e\right)} \epsilon - \frac{1}{2} \, \epsilon e - e^{\left(\frac{-2 \, x}{\epsilon} + 1\right)} + e^{\left(-x + 1\right)}

1-Term Composite Solution: y_{C 0}(x) = -e^{\left(\frac{-2 \, x}{\epsilon} + 1\right)} + e^{\left(-x + 1\right)}
2-Term Composite Solution: y_{C 1}(x) = -\frac{1}{2} \, {\left(x - 1\right)} \epsilon e^{\left(-x + 1\right)} - \frac{1}{2} \, {\left(\frac{4 \, x e^{\left(\frac{-2 \, x}{\epsilon} + 1\right)}}{\epsilon} + e^{\left(\frac{-2 \, x}{\epsilon} + 1\right)} - e\right)} \epsilon - \frac{1}{2} \, \epsilon e - e^{\left(\frac{-2 \, x}{\epsilon} + 1\right)} + e^{\left(-x + 1\right)}

# Exact solution rP = (-1 + sqrt(1 - 2 * epsilon)) / epsilon rM = (-1 - sqrt(1 - 2 * epsilon)) / epsilon yE(x) = (exp(rP * x) - exp(rM * x)) / (exp(rP) - exp(rM)) 

# Interactively plot solution & approximations # Set up controls for interactive plots def showCheckBox(*args, **kwargs): return checkbox(*args, **kwargs) def colorSelector(*args, **kwargs): return color_selector(*args, label = 'Color', hide_box = True, widget = 'jpicker', **kwargs) # The same transparency for each curve transparency = 0.5 lineStyles = ('-', '--', '-.', ':') def lineStyleSelector(*args, **kwargs): return selector(*args, values = lineStyles, label = 'Style', **kwargs) thicknesses = (1, 2, 3) def thicknessSelector(*args, **kwargs): return selector(*args, values = thicknesses, label = 'Thickness', **kwargs) # Number of points on each curve nPoints = int(1001) # Make the plots @interact(layout = [['log10Epsilon'], ['EShown', 'EColor', 'EStyle', 'EThickness'], ['yC0Shown', 'yC0Color', 'yC0Style', 'yC0Thickness'], ['yC1Shown', 'yC1Color', 'yC1Style', 'yC1Thickness']]) def plotSeries(log10Epsilon = slider(-3, log(0.5, 10), default = -2, label = '$\log_{10} \epsilon$'), EShown = showCheckBox(label = 'Exact Solution'), EColor = colorSelector(default = (0, 0, 0)), EStyle = lineStyleSelector(default = '--'), EThickness = thicknessSelector(), yC0Shown = showCheckBox(label = '1-Term Composite Solution'), yC0Color = colorSelector(default = (1, 0, 0)), yC0Style = lineStyleSelector(), yC0Thickness = thicknessSelector(), yC1Shown = showCheckBox(label = '2-Term Composite Solution'), yC1Color = colorSelector(default = (0, 1, 0)), yC1Style = lineStyleSelector(), yC1Thickness = thicknessSelector()): epsilon = 10^log10Epsilon html(r'$$\epsilon = %s$$' % latex(epsilon)) solutionPlots = [] absErrorPlots = [] relErrorPlots = [] # Add boundary points solutionPlots.append(point([(0, 0), (1, 1)], color = 'black', size = 20)) if yC0Shown: solutionPlots.append(plot(yC0(epsilon = epsilon), x, 0, 1, plot_points = nPoints, color = yC0Color, alpha = transparency, linestyle = yC0Style, thickness = yC0Thickness)) absErrorPlots.append(plot((yE - yC0)(epsilon = epsilon), x, 0, 1, plot_points = nPoints, color = yC0Color, alpha = transparency, linestyle = yC0Style, thickness = yC0Thickness)) relErrorPlots.append(plot(((yE - yC0) / yE)(epsilon = epsilon), x, 0, 1, plot_points = nPoints, color = yC0Color, alpha = transparency, linestyle = yC0Style, thickness = yC0Thickness)) if yC1Shown: solutionPlots.append(plot(yC1(epsilon = epsilon), x, 0, 1, plot_points = nPoints, color = yC1Color, alpha = transparency, linestyle = yC1Style, thickness = yC1Thickness)) absErrorPlots.append(plot((yE - yC1)(epsilon = epsilon), x, 0, 1, plot_points = nPoints, color = yC1Color, alpha = transparency, linestyle = yC1Style, thickness = yC1Thickness)) relErrorPlots.append(plot(((yE - yC1) / yE)(epsilon = epsilon), x, 0, 1, plot_points = nPoints, color = yC1Color, alpha = transparency, linestyle = yC1Style, thickness = yC1Thickness)) if EShown: solutionPlots.append(plot(yE(epsilon = epsilon), x, 0, 1, plot_points = nPoints, color = EColor, alpha = transparency, linestyle = EStyle, thickness = EThickness)) show(sum(solutionPlots), axes_labels = (r'$x$', r'$y(x)$')) show(sum(absErrorPlots), axes_labels = (r'$x$', 'Abs. Err.')) show(sum(relErrorPlots), axes_labels = (r'$x$', 'Rel. Err.'))