Figure 10 (page 19):

The function e^{-10x} and 10 collocation points.

Code for Figure 10

Text of the GNU GPL.

main.m

clear('all'); close('all');

n = 10;
[x, A, B, Q] = colloc(n);
a = 10;

nplot = 100;
xplot = linspace(0, 1, nplot)';

z = exp(-a*xplot);

Z = exp(-a*x);

figure(1)
plot(xplot, z,  x, Z, 'o')

table1 = [xplot, z];
table2 = [x, Z];

save expfcn_1.dat table1
save expfcn_2.dat table2

/export/home/jbraw/courses/cbe255/content/util/matlab/colloc.m

function [r, a, b, q] = colloc (n, varargin)

  nargs = nargin;

  if (nargs < 1 || nargs > 3)
    error ('usage: [r, a, b, q] = colloc (n, ''left'', ''right'')');
  end

  if (~ (isnumeric (n) && numel (n) == 1 && round (n) == n))
    error ('colloc: first argument must be an integer scalar');
  end

  if (isnan (n) || isinf (n))
    error ('colloc: NaN is invalid as N');
  end

  if (n < 0)
    error ('colloc: first argument must be non-negative');
  end

  n0 = 0;
  n1 = 0;

  for i = 2:nargs
    s = varargin{i-1};
    if (~ ischar (s))
      error ('colloc: expecting character string argument');
    end

    s = lower (s);
    if (strcmp (s, 'r') || strcmp (s, 'right'))
      n1 = 1;
    elseif (strcmp (s, 'l') || strcmp (s, 'left'))
      n0 = 1;
    else
      error ('colloc: unrecognized argument');
    end
  end

  nt = n + n0 + n1;

  if (nt < 1)
    error ('colloc: the total number of roots must be positive');
  end

  alpha = 0;
  beta = 0;

  %% Compute roots.

  [dif1, dif2, dif3, r] = jcobi (n, n0, n1, alpha, beta);

  %% First derivative weights.

  a = zeros (nt, nt);
  for i = 1:nt
    a(i,:) = dfopr (n, n0, n1, i, 1, dif1, dif2, dif3, r)';
  end

  %% Second derivative weights.

  b = zeros (nt, nt);
  for i = 1:nt
    b(i,:) = dfopr (n, n0, n1, i, 2, dif1, dif2, dif3, r)';
  end

  %% Gaussian quadrature weights.

  id = 3;
  q = dfopr (n, n0, n1, 0, id, dif1, dif2, dif3, r);

end


%% The following routines (JCOBI, DIF, DFOPR, INTRP, AND RADAU)
%% are the same as found in Villadsen, J. and M.L. Michelsen,
%% Solution of Differential Equation Models by Polynomial
%% Approximation, Prentice-Hall (1978) pages 418-420.
%%
%% Cosmetic changes (elimination of arithmetic IF statements, most
%% GO TO statements, and indentation of program blocks) made by:
%%
%% John W. Eaton
%% Department of Chemical Engineering
%% The University of Texas at Austin
%% Austin, Texas 78712
%%
%% June 6, 1987
%%
%% Some error checking additions also made on June 7, 1987
%%
%% Further cosmetic changes made August 20, 1987
%%
%% Translated from Fortran December 14, 2006

function vect = dfopr (n, n0, n1, i, id, dif1, dif2, dif3, root)

%%   Villadsen and Michelsen, pages 133-134, 419
%%
%%   Input parameters:
%%
%%     N      : The degree of the Jacobi polynomial, (i.e. the number
%%              of interior interpolation points)
%%
%%     N0     : Determines whether x = 0 is included as an
%%              interpolation point
%%
%%                n0 = 0  ==>  x = 0 is not included
%%                n0 = 1  ==>  x = 0 is included
%%
%%     N1     : Determines whether x = 1 is included as an
%%              interpolation point
%%
%%                n1 = 0  ==>  x = 1 is not included
%%                n1 = 1  ==>  x = 1 is included
%%
%%     I      : The index of the node for which the weights are to be
%%              calculated
%%
%%     ID     : Indicator
%%
%%                id = 1  ==>  first derivative weights are computed
%%                id = 2  ==>  second derivative weights are computed
%%                id = 3  ==>  Gaussian weights are computed (in this
%%                             case, I is not used).
%%
%%   Output parameters:
%%
%%     DIF1   : vector containing the first derivative
%%              of the node polynomial at the zeros
%%
%%     DIF2   : vector containing the second derivative
%%              of the node polynomial at the zeros
%%
%%     DIF3   : vector containing the third derivative
%%              of the node polynomial at the zeros
%%
%%     VECT   : vector of computed weights

  if (n0 ~= 0 && n0 ~= 1)
    error ('dfopr: n0 not equal to 0 or 1');
  end

  if (n1 ~= 0 && n1 ~= 1)
    error ('dfopr: n1 not equal to 0 or 1');
  end

  if (n < 0)
    error ('dfopr: n less than 0');
  end

  nt = n + n0 + n1;

  if (id ~= 1 && id ~= 2 && id ~= 3)
    error ('dfopr: id not equal to 1, 2, or 3');
  end

  if (id ~= 3)
    if (i < 1)
      error ('dfopr: index less than zero');
    end

    if (i > nt)
      error ('dfopr: index greater than number of interpolation points');
    end
  end

  if (nt < 1)
    error ('dfopr: number of interpolation points less than 1');
  end

%% Evaluate discretization matrices and Gaussian quadrature
%% weights.  Quadrature weights are normalized to sum to one.

  vect = zeros (nt, 1);

  if (id ~= 3)
    for j = 1:nt
      if (j == i)
	if (id == 1)
	  vect(i) = dif2(i)/dif1(i)/2.0;
	else
	  vect(i) = dif3(i)/dif1(i)/3.0;
	end
      else
	y = root(i) - root(j);
	vect(j) = dif1(i)/dif1(j)/y;
	if (id == 2)
	  vect(j) = vect(j)*(dif2(i)/dif1(i) - 2.0/y);
	end
      end
    end
  else
    y = 0.0;

    for j = 1:nt

      x = root(j);
      ax = x*(1.0 - x);

      if (n0 == 0)
	ax = ax/x/x;
      end

      if (n1 == 0)
	ax = ax/(1.0 - x)/(1.0 - x);
      end

      vect(j) = ax/dif1(j)^2;
      y = y + vect(j);

    end

    vect = vect/y;

  end

end

function [dif1, dif2, dif3, root] = jcobi (n, n0, n1, alpha, beta)

%%   Villadsen and Michelsen, pages 131-132, 418
%%
%%   This subroutine computes the zeros of the Jacobi polynomial
%%
%%      (ALPHA,BETA)
%%     P  (X)
%%      N
%%
%%   Input parameters:
%%
%%     N      : The degree of the Jacobi polynomial, (i.e. the number
%%              of interior interpolation points)
%%
%%     N0     : Determines whether x = 0 is included as an
%%              interpolation point
%%
%%                N0 = 0  ==>  x = 0 is not included
%%                N0 = 1  ==>  x = 0 is included
%%
%%     N1     : Determines whether x = 1 is included as an
%%              interpolation point
%%
%%                N1 = 0  ==>  x = 1 is not included
%%                N1 = 1  ==>  x = 1 is included
%%
%%     ALPHA  : The value of alpha in the description of the jacobi
%%              polynomial
%%
%%     BETA   : The value of beta in the description of the Jacobi
%%              polynomial
%%
%%     For a more complete explanation of alpha and beta, see Villadsen
%%     and Michelsen, pages 57 to 59
%%
%%   Output parameters:
%%
%%     ROOT   : vector containing the n + n0 + n1 zeros of the node
%%              polynomial used in the interpolation routine
%%
%%     DIF1   : vector containing the first derivative
%%              of the node polynomial at the zeros
%%
%%     DIF2   : vector containing the second derivative
%%              of the node polynomial at the zeros
%%
%%     DIF3   : vector containing the third derivative
%%              of the node polynomial at the zeros

  if (n0 ~= 0 && n0 ~= 1)
    error ('jcobi: n0 not equal to 0 or 1');
  end

  if (n1 ~= 0 && n1 ~= 1)
    error ('jcobi: n1 not equal to 0 or 1');
  end

  if (n < 0)
    error ('jcobi: n less than 0');
  end

  nt = n + n0 + n1;

  if (nt < 1)
    error ('jcobi: number of interpolation points less than 1');
  end

  dif1 = zeros (nt, 1);
  dif2 = zeros (nt, 1);
  dif3 = zeros (nt, 1);
  root = zeros (nt, 1);

%% First evaluation of coefficients in recursion formulas.
%% recursion coefficients are stored in dif1 and dif2.

  ab = alpha + beta;
  ad = beta - alpha;
  ap = beta*alpha;
  dif1(1) = (ad/(ab + 2.0) + 1.0)/2.0;
  dif2(1) = 0.0;

  if (n >= 2)
    for i = 2:n

      z1 = i - 1.0;
      z = ab + 2*z1;
      dif1(i) = (ab*ad/z/(z + 2.0) + 1.0)/2.0;

      if (i == 2)
	dif2(i) = (ab + ap + z1)/z/z/(z + 1.0);
      else
	z = z*z;
	y = z1*(ab + z1);
	y = y*(ap + y);
	dif2(i) = y/z/(z - 1.0);
      end

    end
  end

%% Root determination by newton method with suppression of
%% previously determined roots.

  x = 0.0;

  for i = 1:n

    done = false;

    while (~ done)

      xd = 0.0;
      xn = 1.0;
      xd1 = 0.0;
      xn1 = 0.0;

      for j = 1:n
	xp = (dif1(j) - x)*xn  - dif2(j)*xd;
	xp1 = (dif1(j) - x)*xn1 - dif2(j)*xd1 - xn;
	xd = xn;
	xd1 = xn1;
	xn = xp;
	xn1 = xp1;
      end

      zc = 1.0;
      z = xn/xn1;

      if (i ~= 1)
	for j = 2:i
	  zc = zc - z/(x - root(j-1));
	end
      end

      z = z/zc;
      x = x - z;

      if (abs(z) <= 1.0e-09)
	done = true;
      end

    end

    root(i) = x;
    x = x + 0.0001;

  end

%% Add interpolation points at x = 0 and/or x = 1.

  nt = n + n0 + n1;

  if (n0 ~= 0)
    for i = 1:n
      j = n + 1 - i;
      root(j+1) = root(j);
    end
    root(1) = 0.0;
  end

  if (n1 == 1)
    root(nt) = 1.0;
  end


%% Use recursion formulas to evaluate derivatives of node polynomial
%%
%%                   N0     (ALPHA,BETA)           N1
%%     P  (X)  =  (X)   *  P (X)         *  (1 - X)
%%      NT                   N
%%
%% at the interpolation points.

  for i = 1:nt
    x = root(i);
    dif1(i) = 1.0;
    dif2(i) = 0.0;
    dif3(i) = 0.0;
    for j = 1:nt
      if (j ~= i)
	y = x - root(j);
	dif3(i) = y*dif3(i) + 3.0*dif2(i);
	dif2(i) = y*dif2(i) + 2.0*dif1(i);
	dif1(i) = y*dif1(i);
      end
    end
  end

end

Data files

expfcn_1.dat

# Created by Octave 3.8.1, Thu Aug 21 10:29:16 2014 CDT 
# name: table1
# type: matrix
# rows: 100
# columns: 2
 0 1
 0.0101010101010101 0.9039239022952824
 0.0202020202020202 0.8170784211407313
 0.0303030303030303 0.7385767149187981
 0.04040404040404041 0.6676171462938303
 0.05050505050505051 0.6034750961171595
 0.06060606060606061 0.5454955638202434
 0.07070707070707072 0.4930864787331597
 0.08080808080808081 0.4457126540255176
 0.09090909090909091 0.402890321529133
 0.101010101010101 0.3641821916336149
 0.1111111111111111 0.3291929878079056
 0.1212121212121212 0.2975654101475653
 0.1313131313131313 0.2689764867286835
 0.1414141414141414 0.2431342755094668
 0.1515151515151515 0.2197748831002536
 0.1616161616161616 0.1986597699584707
 0.1717171717171717 0.179573314489944
 0.1818181818181818 0.1623206111818482
 0.1919191919191919 0.1467254802824514
 0.202020202020202 0.132628668703063
 0.2121212121212121 0.1198862237703009
 0.2222222222222222 0.1083680232218959
 0.2323232323232323 0.09795644643476188
 0.2424242424242424 0.08854517331628878
 0.2525252525252525 0.08003809859347186
 0.2626262626262627 0.07234835041290566
 0.2727272727272728 0.06539740322986018
 0.2828282828282829 0.0591142759275133
 0.2929292929292929 0.05343480697775794
 0.303030303030303 0.04830099924173013
 0.3131313131313131 0.04366042771934617
 0.3232323232323233 0.0394657041999525
 0.3333333333333334 0.03567399334725237
 0.3434343434343435 0.03224657527690432
 0.3535353535353536 0.02914845015995792
 0.3636363636363636 0.02634798081444872
 0.3737373737373738 0.02381656963539773
 0.3838383838383839 0.02152836656411605
 0.393939393939394 0.01946000511467905
 0.4040404040404041 0.01759036376194685
 0.4141414141414142 0.01590035025449252
 0.4242424242424243 0.01437270664990267
 0.4343434343434344 0.01299183308152538
 0.4444444444444445 0.01174362845702136
 0.4545454545454546 0.01061534646197667
 0.4646464646464647 0.009595465398126373
 0.4747474747474748 0.008673570527013745
 0.4848484848484849 0.007840247717611619
 0.494949494949495 0.007086987311865175
 0.5050505050505051 0.006406097226458321
 0.5151515151515152 0.00579062440342319
 0.5252525252525253 0.005234283807468586
 0.5353535353535354 0.004731394244968012
 0.5454545454545455 0.004276820349208926
 0.5555555555555556 0.003865920139472808
 0.5656565656565657 0.003494497618434179
 0.5757575757575758 0.003158759923816596
 0.5858585858585859 0.002855278596750249
 0.595959595959596 0.00258095457131468
 0.6060606060606061 0.002332986527749615
 0.6161616161616162 0.002108842286165751
 0.6262626262626263 0.001906232948636252
 0.6363636363636365 0.001723089525615121
 0.6464646464646465 0.001557541807998149
 0.6565656565656566 0.001407899269073737
 0.6666666666666667 0.001272633801339807
 0.6767676767676768 0.001150364111899958
 0.686868686868687 0.001039841617089057
 0.696969696969697 0.000939937692288177
 0.7070707070707072 0.0008496321467275512
 0.7171717171717172 0.0007680028055854866
 0.7272727272727273 0.0006942160929985581
 0.7373737373737375 0.0006275185198194412
 0.7474747474747475 0.0005672289891977492
 0.7575757575757577 0.0005127318414106374
 0.7676767676767677 0.0004634705669189497
 0.7777777777777778 0.0004189421234483841
 0.787878787878788 0.0003786917990633349
 0.797979797979798 0.0003423085687765509
 0.8080808080808082 0.0003094208972776129
 0.8181818181818182 0.0002796929449188878
 0.8282828282828284 0.0002528211382155403
 0.8383838383838385 0.0002285310698385263
 0.8484848484848485 0.0002065746964441565
 0.8585858585858587 0.0001867278057252652
 0.8686868686868687 0.0001687877268182172
 0.8787878787878789 0.0001525712606850728
 0.888888888888889 0.0001379128093365619
 0.8989898989898991 0.0001246626847920102
 0.9090909090909092 0.0001126855805078007
 0.9191919191919192 0.0001018591896650205
 0.9292929292929294 9.20729562066405e-05
 0.9393939393939394 8.322694587016918e-05
 0.9494949494949496 7.523082568708149e-05
 0.9595959595959597 6.800294152796293e-05
 0.9696969696969697 6.14694842735142e-05
 0.9797979797979799 5.556373609659339e-05
 0.9898989898989899 5.022538915853798e-05
 1 4.539992976248485e-05


expfcn_2.dat

# Created by Octave 3.8.1, Thu Aug 21 10:29:16 2014 CDT 
# name: table2
# type: matrix
# rows: 10
# columns: 2
 0.01304673574141415 0.8776851423937573
 0.06746831665550772 0.5093177639484284
 0.1602952158504878 0.201301366395121
 0.2833023029353764 0.05883472545961437
 0.4255628305091844 0.01418417595529056
 0.5744371694908156 0.003200744964359464
 0.7166976970646236 0.0007716519352784012
 0.8397047841495122 0.0002255321490136947
 0.9325316833444922 8.913871256036119e-05
 0.9869532642585859 5.172689791542233e-05