forked from MrFelixP/PoCET
-
Notifications
You must be signed in to change notification settings - Fork 1
/
PoCETquadRules.m
52 lines (46 loc) · 1.7 KB
/
PoCETquadRules.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
function quadrules = PoCETquadRules(PoCETsys)
% q = PoCETquadRules(PoCETsys)
% returns a struct containing the roots, weights, normalized weights,
% polynomials, and normalized polynomials required for integration via
% gaussian quadrature. All polynomials are already evaluated at the roots.
opts = PoCETsys.pce.options;
state_names = {PoCETsys.states(:).name};
if ~isempty(PoCETsys.parameters(1).name)
param_names = {PoCETsys.parameters(:).name};
else
param_names = {};
end
qp = ceil(opts.hpo*opts.order/2 + 0.5);
% preallocation
rj = zeros(opts.n_xi,qp);
wj = zeros(opts.n_xi,qp);
PHI_rj = zeros(opts.order+1,qp,opts.n_xi);
% calculate gaussian quadrature roots & weights and evaluate polynomials
for i = 1:opts.n_xi
statenum = find(ismember(state_names,PoCETsys.pce.vars{i}));
paramnum = find(ismember(param_names,PoCETsys.pce.vars{i}));
if ~isempty(statenum)
curvar = PoCETsys.states(statenum);
else
curvar = PoCETsys.parameters(paramnum);
end
if strcmp(curvar.dist,'normal')
[rj(i,1:qp), wj(i,1:qp)] = gauss_hermite_quadrature(qp);
wj(i,:) = wj(i,:)/sqrt(2*pi);
PHI_rj(:,:,i) = calc_hermite_polynomials(rj(i,:),opts.order,qp);
elseif strcmp(curvar.dist,'uniform')
[rj(i,1:qp), wj(i,1:qp)] = gauss_legendre_quadrature(qp);
wj(i,:) = wj(i,:)/2;
PHI_rj(:,:,i) = calc_legendre_polynomials(rj(i,:),opts.order,qp);
elseif strcmp(curvar.dist,'beta') || strcmp(curvar.dist,'beta4')
alf = curvar.data(1);
bet = curvar.data(2);
[rj(i,1:qp), wj(i,1:qp)] = gauss_jacobi_quadrature(qp,bet-1,alf-1);
wj(i,:) = wj(i,:)/(2^(alf+bet-1)*beta(alf,bet));
PHI_rj(:,:,i) = calc_jacobi_polynomials(rj(i,:),bet-1,alf-1,opts.order,qp);
end
end
quadrules.rj = rj;
quadrules.wj = wj;
quadrules.PHI_rj = PHI_rj;
end