echo off
%DIFEQDEM singularly perturbed ODE solved by collocation at Gauss points

%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc. 
%   latest change:     June 3, '90
%   cb : 9 may '95 (use .' instead of '); 28sep96; 8aug97
%   cb : 5 oct 97 (replace use of ANS, to help compilation)
%   $Revision: 1.10 $

home, close, echo on
% We consider the nonlinear singularly perturbed problem
%
%  epsilon D^2g(x) + (g(x))^2 = 1  on  [0..1]
%                   Dg(0) = g(1) = 0 .
%

   pause                % Touch any key to continue

% We seek an approximate solution by collocation from C^1 piecewise
% cubics with a suitable break sequence; e.g.,

breaks = [0:4]/4;

% gives the corresponding knot sequence

knots = augknt(breaks,4,2) 

   pause                % Touch any key to continue


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

k = 4;

% This implies that we have

n = length(knots)-k

% degrees of freedom. This fits nicely with the fact that we expect to
% collocate at two points per polynomial piece, for a total of 8 conditions,
% bringing us to 10 conditions altogether once we add the two side conditions.

   pause                % Touch any key to continue

%    We choose two Gauss points for each interval. For the `standard'
% interval [-1 .. 1]/2 of unit length, these are the two points

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

%    From this, we obtain the whole collection of collocation points by

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

   pause                % Touch any key to continue

%    The numerical problem we want to solve is to find a pp  y  of the given
% order and with the given knots that satisfies the (nonlinear) system
%
%                           Dy(0)  =  0
%      (y(x))^2 + epsilon D^2y(x)  =  1   for  x in  colpnts
%                            y(1)  =  0

% If  y  is our current approximation to the solution, then the linear
% problem for the better (?) solution  z  by Newton's method reads
%
%                           Dz(0)  =  0
%    w_0(x)z(x) + epsilon D^2z(x)  =  b(x)   for  x in  colpnts
%                            z(1)  =  0
%
% with  w_0(x) := 2y(x),  b(x) := (y(x))^2 + 1 .

   pause                % Touch any key to continue

% In fact, by choosing   w_0(1) := 1, w_1(0) := 1 , and
%
%     w_2(x) := epsilon,    w_1(x) := 0    for  x  in  colpnts
%
% and choosing all other values of  w_0, w_1, w_2, b  not yet specified
% to be zero, we can give our system the uniform shape

%  w_0(x)z(x) + w_1(x)Dz(x) + w_2(x)D^2z(x) = b(x)   for  x  in  points

% with

points = [0,colpnts,1];

   pause                % Touch any key to continue

%   This system converts into one for the B-spline coefficients of its
% solution  z . For this, we need the zeroth, first, and second derivative
% at every  x  in  points  and for every relevant B-spline. These values
% are supplied by the M-file  spcol .

% Here is the essential part of the online help for spcol:

%SPCOL B-spline collocation matrix.
%
%   SPCOL(KNOTS,K,TAU)  returns the matrix   (D^m(i)B_j(TAU_i)) , with B_j
%   the j-th B-spline of order K for the knot sequence KNOTS, TAU a
%   sequence of points, both assumed to be  n o n d e c r e a s i n g ,
%   and   m(i) := #{ j<i : TAU(j) = TAU(i) }.
%

% We use  spcol  to supply the matrix
colmat = spcol(knots,k,sort([points points points]));

% From this, we get the collocation matrix by combining the row triple
% associated with  x  using the weights  w_0(x), w_1(x), w_2(x)  to get the
% row corresponding to  x  of the matrix of our linear system.

   pause                % Touch any key to continue


% For this, we need to choose  epsilon . The problem is already very
% difficult for  epsilon = .001 , so we choose the mild

epsilon = .1;

%     We also need a current approximation  y . Initially, we get it
% by interpolating some reasonable initial guess from our pp space at all
% the  points .  We use the function
%   ()^2 - 1
% which does satisfy the end conditions, 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 colpnts.*colpnts-1 0].';
y = spmak(knots,coefs.');
% We plot the result (it should be exactly  ()^2-1 ), to be sure:
fnplt(y,'g'), grid off, axis(axis)
title('initial guess is ()^2-1'), pause
hold on
% That's ok, then.

   pause                % Touch any key to continue

%    We can now complete the construction and solution of the linear
% system for the improved approximate solution  z  from our current
% guess  y . In fact, with the initial guess  y  available, we now set
% up an iteration, to be terminated when the change  z-y  is `small
% enough' .  We choose the tolerance

tolerance = 6.e-9;

   pause                % Touch any key to continue

difeqite

% that looks like quadratic convergence, as expected.

   pause                % Touch any key to continue

   difeqset

%  If we now decrease  epsilon , we create more of a boundary layer near
% the right endpoint, and this calls for a nonuniform mesh. We use  newknt
% to construct an appropriate (finer) mesh from the current approximation,
% then divide  epsilon  by 3 and repeat the above calculation.

while 1
   echo on
   difeqset  % generates new breaks, new collocation points, and
             % constructs a new approximate solution by interpolating the
             % current approximation from the new pp space.

   % We plot the result (it should be close to our current solution) to be sure:
   fnplt(y,'g')
   title('new initial guess for new value of epsilon'), xlabel(''), pause
   epsilon = epsilon/3
   difeqite  % runs the Newton iteration

   echo off
   fprintf(' If you want to decrease  epsilon  further,\n')
   fprintf(' type any digit, then RETURN. ... (then wait) ...')
     answer = input(' Otherwise, just hit RETURN :  ');
   if isempty(answer), break,end
end
echo on

% Here is the final distribution of breaks

breaks = fnbrk(fn2fm(z,'pp'),'b');
for j = 1:length(breaks)
   plot(breaks(j)*[1 1],[0,-1],'r')
end
xlabel('... and iterates, and final distribution of breaks.'), pause

hold off, echo off, close
