-
Notifications
You must be signed in to change notification settings - Fork 2
/
Smolyak_Polynomial.m
133 lines (116 loc) · 7.62 KB
/
Smolyak_Polynomial.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
% Smolyak_Polynomial.m is a routine that constructs the multidimensional
% basis functions of Smolyak polynomial of the approximation level that
% corresponds to the previously constructed Smolyak (multidimensional) grid
% points; see "Smolyak method for solving dynamic economic models: Lagrange
% interpolation, anisotropic grid and adaptive domain" by Kenneth L. Judd,
% Lilia Maliar, Serguei Maliar and Rafael, (2014). Journal of Economic
% Dynamics and Control 44, 92–123 (henceforth, JMMV (2014)).
%
% This version: Novenber 5, 2014. First version: May 30, 2011.
% -------------------------------------------------------------------------
% Inputs: "points" is the matrix of points in which the polynomial basis
% functions must be evaluated; numb_pts-by-d
% "d" is the number of dimensions (state variables)
% "Smol_elem" is the matrix of the subindices of the Smolyak
% unidimensional elements; these elements can be either
% isotropic (produced by Smolyak_Elem_Isotrop.m) or
% anisotropic (produced by Smolyak_Elem_Anisotrop.m);
% in the former case, the indices i1,...,id that jointly
% satisfy the Smolyak rule, d<=|i|<=d+mu, where
% |i|=i1+i2+...+id; see JMMV (2014), Section 3.2.3;
% in the later case, they are a subset of the above
% indices
%
% Output: "Smol_bases" is the matrix of multidimensional basis functions of
% Smolyak polynomial of the given level of approximation,
% evaluated in data matrix "points"
% -------------------------------------------------------------------------
% Copyright © 2014 by Lilia Maliar, Serguei Maliar and Rafael Valero. All
% rights reserved. The code may be used, modified and redistributed under
% the terms provided in the file "License_Agreement.txt".
% -------------------------------------------------------------------------
function Smol_bases = Smolyak_Polynomial(points,d,mu,Smol_elem)
% Smolyak polynomial is given by the sum of multidimensional basis functions,
% multiplied by the coefficients; see formula (15) in JMMV (2014). By
% convention, the first basis function is given by 1 (unity).
% Unidimensional basis functions are given by Chebyshev polynomial bases;
% in JMMV (2014), a unidimensional Chebyshev polynomial basis function of
% degree n-1 is denoted by "phi_n", i.e., has a subindex n and we follow
% this notation here
% 1. Construct the unidimensional basis functions and evaluate them in
% all the points of matrix "points"
% -------------------------------------------------------------------------
i_max = mu+1; % The maximum subindex of unidimensional
% set S_i whose points are used to
% construct Smolyak grid; e.g., for mu=1,
% we consider up to S_i_max={0,-1,1} where
% i_max=mu+1=1+1=2
% Compute the number of elements in the i_max-th unidimensional set of
% elements S_i_max; this coincides with a maximum subindex of elements
% (unidimensional grid point or unidimensional basis function);
% see Section 2.2.1 in JMMV (2014)
m_i_max(i_max==1) = 1; % If i_max=1, then m(i_max)=1, i.e., set
% S_1={0} and the maximum subindex is 1
m_i_max(i_max>1) = 2.^(i_max(i_max>1)-1) + 1;
% If i_max>1, then m(i_max)= 2^(i_max-1)+1;
% e.g., for S_2={0,-1,1}, the maximum
% subindex is 3
numb_pts = size(points,1); % Compute the number of points (rows),
% "numb_pts" in the matrix of points
% "points", in which the polynomial bases
% must be evaluated
phi = ones(numb_pts,d,m_i_max);
% Allocate memory to a matrix of the
% unidimensional polynomial bases "phi_n",
% evaluated in all the points
% For a polynomial bases "phi_n" with n=1, we have phi_n(x)=1 for all x;
% our phi(:,:,1) is a matrix of ones of size numb_pts-by-d by the above
% construction
phi(:,:,2) = points; % For a polynomial bases "phi_n" with n=2,
% we have phi_n(x) is x; evaluating it in
% all the points gives us matrix "points";
% numb_pts-by-d
for j = 3:m_i_max % For polynomial bases "phi_n", from n=3 to
% n=m_i_max, ...
phi(:,:,j) = 2.* phi(:,:,2).*phi(:,:,j-1) - phi(:,:,j-2);
% Use the recurrence formula to compute the
% Chebyshev polynomial basis functions of
% the degrees from 2 to m_i_max-1
end
% 2. Form the multidimensional polynomial bases of Smolyak polynomial of the
% required level of Smolyak approximation; see JMMV (2014), Sections 3.3.3
% and 3.4.2 for examples
% ----------------------------------------------------------------------
Smol_bases = []; % Initially, the matrix of multidimensional
% polynomial bases is empty
numb_terms = size(Smol_elem,1);
% Compute the number of terms (i.e., multi-
% dimensional polynomial bases) in Smolyak
% polynomial
for jt = 1:numb_terms % For each term of Smolyak polynomial, ...
index_row = Smol_elem(jt,:);
% Identify the subindices of the unidimensional
% basis function that constitute a jt-th multi-
% dimensional basis function; this is a jt-th
% row of matrix "Smol_elem"; 1-by-d
product = ones(numb_pts,1);
% Initialize a vector of products of unidimen-
% sional basis functions; numb_pts-by-1
for jd = 1:d % For each dimension (state variable), ...
n = index_row(jd); % A subindex of a unidimensional basis function
% phi_n in a dimension jd is denoted n
if n ~= 1; % If the subindex of unidimensional basis
% function is not equal to unity, ...
product = product.*phi(:,jd,n);
% Compute the product of basis functions
% Otherwise, i.e., if n = 1, there is no need to compute the
% product of unidimensional basis functions, as it's equal to unity
end
end
Smol_bases = cat(2,Smol_bases,product);
% Attach to the previously obtained matrix of
% multidimensional basis functions a new
% product of unidimensional basis functions;
% e.g., for mu=1 and d=2, basis_bs is of
% size numb_pts-by-5
end