% DIFEQDEM	Demonstrate solution of ODE by collocation.
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
%                                       latest change:     June 3, '90
clg;clc;home;echo on
% We consider the nonlinear singularly perturbed problem
%
%  epsilon g''(x) + (g(x))^2 = 1  on  [0,1]
%                   g'(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 breakpoint sequence; e.g.,

breaks = [0:4]/4;

% hence with the knot sequence

knots = sort([0 0 breaks breaks 1 1])

% which has the same effect as the statement 
%    knots = augknt(breaks,4,2)
   pause                % Touch any key to continue
echo off       % for a quick check
zcheck(knots-augknt(breaks,4,2));echo on


%     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;
(breaks(2:ninterv+1)+breaks(1:ninterv))/2;
ans'*[1 1] + diff(breaks)'*gauss;
colpnts = sort(ans(:)');

   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 which satisfies the (nonlinear) system
%     
%                           y'(0)  =  0
%       (y(x))^2 + epsilon y''(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
%     
%                           z'(0)  =  0
%     w_0(x)z(x) + epsilon z''(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)z'(x) + w_2(x)z''(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 .

help spcol

   pause                % Touch any key to continue

% 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);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

echo off %  (for a check of newknt etc. that will take a moment)
typ =[-0.97210037583603;
      -0.97138283316330;
      -0.96077877735438;
      -0.94305590383320;
      -0.89729018696592;
      -0.83892108187956;
      -0.70191870276411;
      -0.53656329907016;
      -0.17862024561499;
       0];
zcheck(typ'-fnval(y,points));

   difeqset

echo off       % for a quick check
   typ=[-0.97209501816518;
        -0.97162227476240;
        -0.96470968097732;
        -0.95393233409872;
        -0.93022030933061;
        -0.90278210186314;
        -0.84585242261528;
        -0.78262837036445;
        -0.64997520321340;
        -0.49555677657819;
        -0.16116850509888;
                        0];
   zcheck(fnval(y,points)'-typ);
echo on

%  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,'x');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 = ppbrk(sp2pp(z));
for j = 1:length(breaks)
   plot(breaks(j)*[1 1],[0,-1])
end

hold off
