Perturbation Theory Problems with bvp4c

- 2 minutes read - 300 words

I have been watching Nathan Kutz' lectures on Coursera. One change he made to the course since I took AMATH 581 at University of Washington was introducing the MATLAB function bvp4c. I immediately realized that this would be nice for solving boundary layer problems that arise in asymptotics.

Following my life philosophy of doing the dumbest thing that could possibly work, I tried implementing Nathan’s code for a single-layer boundary layer problem from Holmes, Chapter 2:

$$ \begin{align} \text{DE:} \qquad & \epsilon y'' + 2y' + 2y = 0\ \text{BC:} \qquad & y(0)=0, \quad y(1)=1 \end{align} $$

The code for this DE looks like:

clear all; close all;clc
epsilonVec = [1,1e-1,1e-2];
BS={};
for i=1:length(epsilonVec);
	epsilon = epsilonVec(i);
	singular_rhs_anon=@(x,y) [y(2); (-2*y(2)-2*y(1))/epsilon];
	singular_bc_anon=@(yl,yr) [yl(1)-0; yr(1)-1];

	init = bvpinit(linspace(0,1,20), [0,0]);
	sol = bvp4c(singular_rhs_anon, ...
	            singular_bc_anon, init);
	x = linspace(0,1,1000);
	BS{i}=deval(sol, x);
end
plot(x,BS{1}(1,:),x,BS{2}(1,:),x,BS{3}(1,:));
legend('\epsilon=1','\epsilon=1e-1','\epsilon=1e-2');

Here, the rhs function works exactly like a standard ODE45 call(turn the second order DE into a system of first order DE; coefficient on $y''$ term should be 1). The bc function specifies the right and left bounds. The guess given to bvpinit is critical, but it seems that $(0,0)$ worked okay here. Also note that the granularity might need to be turned up for more difficult problems. Finally, note the use of anonymous functions to define a function without a separate MATLAB .m file.

This generates plots like the following: Boundary Layer with BVP4c

This is a pretty slick way to solve boundary layer problems!

There’s a lot of existing stuff out there on this. A quick google search for ‘perturbation problems with bvp4c’ turned up: