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