Solving A Nonlinear ODE

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”.

Problem

Consider the nonlinear singularly perturbed problem:

$\begin{matrix} ε D^{2} g (x) + {(g (x))}^{2} = 1 & o n & [0..1] \end{matrix}$

$D g (0) = g (1) = 0$

Approximation Space

Seek an approximate solution by collocation from C¹ 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;

i.e., 10 degrees of freedom.

Discretization

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(:).';

Numerical Problem

With this, the numerical problem you want to solve is to find $y \in S_{4, k n o t s}$ that satisfies the nonlinear system

$\begin{matrix} D y (0) = 0 \\ (^{y (x)) 2} + ε D^{2} y (x) = 1 for x \in colsites \\ y (1) = 0 \end{matrix}$

Linearization

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{matrix} D z (0) = 0 \\ w_{0} (x) z (x) + ε D^{2} z (x) = b (x) for x \in colsites \\ z (1)=0 \end{matrix}$

with w₀(x)=2y(x),b(x)=(y(x))²+1. In fact, by choosing

$\begin{array}{l} w_{0} (1) : = 1, w_{1} (0) : = 1 \\ w_{1} (x) : = 0, w_{2} (x) : = ε for x \in colsites \end{array}$

and choosing all other values of w₀,w₁,w₂, b not yet specified to be zero, you can give your system the uniform shape

$\begin{matrix} w_{0} (x) z (x) + w_{1} (x) D z (x) + w_{2} (x) D^{2} z (x) = b (x), & for & x \in sites \end{matrix}$

with

sites = [0,colsites,1];

Linear System to Be Solved

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₀(x),w₁(x),w₂(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²–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

Iteration

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.