This section discusses these aspects of a nonlinear ODE problem:

You can run this example: “Solving a Nonlinear ODE with a Boundary Layer by Collocation”.

Consider the nonlinear singularly perturbed problem:

$$\begin{array}{ccc}\epsilon {D}^{2}g\left(x\right)+{\left(g\left(x\right)\right)}^{2}=1& on& \left[\mathrm{0..1}\right]\end{array}$$

$$Dg\left(0\right)=g\left(1\right)=0$$

Seek an approximate solution by collocation from *C*^{1}^{ } piecewise
cubics with a suitable break sequence; for instance,

breaks = (0:4)/4;

Because cubics are of order 4, you have

k = 4;

Obtain the corresponding knot sequence as

knots = augknt(breaks,k,2);

This gives a quadruple knot at both 0 and 1, which is consistent with the fact that you have cubics, i.e., have order 4.

This implies that you have

n = length(knots)-k; n = 10;

You collocate at two sites per polynomial piece, i.e., at eight sites altogether. This, together with the two side conditions, gives us 10 conditions, which matches the 10 degrees of freedom.

Choose the two Gaussian sites for each interval. For the *standard* interval
[–0.5,0.5] of length 1, these are the two sites

gauss = .5773502692*[-1/2; 1/2];

From this, you obtain the whole collection of collocation sites by

ninterv = length(breaks)-1; temp = ((breaks(2:ninterv+1)+breaks(1:ninterv))/2); temp = temp([1 1],:) + gauss*diff(breaks); colsites = temp(:).';

With this, the numerical problem you want to solve is to find $$y\in {S}_{4,knots}$$ that satisfies the nonlinear system

$$\begin{array}{c}Dy(0)=0\\ {(}^{y}+\epsilon {D}^{2}y(x)=1\text{for}x\text{}\in \text{colsites}\\ y(1)=0\end{array}$$

If *y* is your current approximation to the
solution, then the linear problem for the supposedly better solution *z* by Newton's method reads

$$\begin{array}{c}Dz(0)=0\\ {w}_{0}(x)z(x)+\epsilon {D}^{2}z(x)=b(x)\text{for}x\text{}\in \text{colsites}\\ z\text{(1)=0}\end{array}$$

with *w*_{0}(*x*)=2*y*(*x*),*b*(*x*)=(*y*(*x*))^{2}+1.
In fact, by choosing

$$\begin{array}{l}{w}_{0}(1):=1,\text{}{w}_{1}(0):=1\\ {w}_{1}(x):=0,\text{}{w}_{2}(x):=\epsilon \text{for}x\in \text{colsites}\end{array}$$

and choosing all other values of *w*_{0},*w*_{1},*w*_{2}, *b* not
yet specified to be zero, you can give your system the uniform shape

$$\begin{array}{ccc}{w}_{0}\left(x\right)z\left(x\right)+{w}_{1}\left(x\right)Dz\left(x\right)+{w}_{2}\left(x\right){D}^{2}z\left(x\right)=b\left(x\right),& \text{for}& x\text{}\in \text{sites}\end{array}$$

with

sites = [0,colsites,1];

Because *z*∊*S*_{4,knots}, convert this
last system into a system for the B-spline coefficients of *z*. This
requires the values, first, and second derivatives at every
*x*∊*sites* and for all the relevant B-splines. The
command `spcol`

was expressly written for this purpose.

Use `spcol`

to supply the matrix

colmat = ... spcol(knots,k,brk2knt(sites,3));

From this, you get the collocation matrix by combining the row
triple of `colmat`

for *x* using
the weights *w*_{0}(*x*),*w*_{1}(*x*),*w*_{2}(*x*)
to get the row for *x* of the actual matrix. For
this, you need a current approximation *y*. Initially,
you get it by interpolating some reasonable initial guess from your
piecewise-polynomial space at the `sites`

. Use the
parabola *x*^{2}–1,
which satisfies the end conditions as the initial guess, and pick
the matrix from the full matrix `colmat`

. Here it
is, in several cautious steps:

intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3],:); coefs = intmat\[0 colsites.*colsites-1 0].'; y = spmak(knots,coefs.');

Plot the initial guess, and turn hold on for subsequent plotting:

fnplt(y,'g'); legend('Initial Guess (x^2-1)','location','NW'); axis([-0.01 1.01 -1.01 0.01]); hold on

You can now complete the construction and solution of the linear
system for the improved approximate solution *z* from
your current guess *y*. In fact, with the initial
guess *y* available, you now set up an iteration,
to be terminated when the change *z*–*y* is *small
enough*. Choose a relatively mild ε
= .1.

tolerance = 6.e-9; epsilon = .1; while 1 vtau = fnval(y,colsites); weights=[0 1 0; [2*vtau.' zeros(n-2,1) repmat(epsilon,n-2,1)]; 1 0 0]; colloc = zeros(n,n); for j=1:n colloc(j,:) = weights(j,:)*colmat(3*(j-1)+(1:3),:); end coefs = colloc\[0 vtau.*vtau+1 0].'; z = spmak(knots,coefs.'); fnplt(z,'k'); maxdif = max(max(abs(z.coefs-y.coefs))); fprintf('maxdif = %g\n',maxdif) if (maxdif<tolerance), break, end % now reiterate y = z; end legend({'Initial Guess (x^2-1)' 'Iterates'},'location','NW');

The resulting printout of the errors is:

maxdif = 0.206695 maxdif = 0.01207 maxdif = 3.95151e-005 maxdif = 4.43216e-010

If you now decrease ε, you create more of a boundary layer near the right endpoint, and this calls for a nonuniform mesh.

Use `newknt`

to construct
an appropriate finer mesh from the current approximation:

knots = newknt(z, ninterv+1); breaks = knt2brk(knots); knots = augknt(breaks,4,2); n = length(knots)-k;

From the new break sequence, you generate the new collocation site sequence:

ninterv = length(breaks)-1; temp = ((breaks(2:ninterv+1)+breaks(1:ninterv))/2); temp = temp([1 1], :) + gauss*diff(breaks); colpnts = temp(:).'; sites = [0,colpnts,1];

Use `spcol`

to supply the
matrix

colmat = spcol(knots,k,sort([sites sites sites]));

and use your current approximate solution `z`

as
the initial guess:

intmat = colmat([2 1+(1:(n-2))*3,1+(n-1)*3],:); y = spmak(knots,[0 fnval(z,colpnts) 0]/intmat.');

Thus set up, divide ε by 3 and repeat the earlier calculation, starting with the statements

tolerance=1.e-9; while 1 vtau=fnval(y,colpnts); . . .

Repeated passes through this process generate a sequence of solutions, for ε = 1/10, 1/30, 1/90, 1/270, 1/810. The resulting solutions, ever flatter at 0 and ever steeper at 1, are shown in the example plot. The plot also shows the final break sequence, as a sequence of vertical bars. To view the plots, run the example “Solving a Nonlinear ODE with a Boundary Layer by Collocation”.

In this example, at least, `newknt`

has performed
satisfactorily.