Skip to content

Commit

Permalink
updates docs and cleans up files for differentiation and integration
Browse files Browse the repository at this point in the history
of Legendre series. this fixes the behavior of intpol. Previously,
intpol set the constant in a way that assumed you were resampling on
the legendre nodes and cancelled the highest degree term. Now, this behvavior
is only used for generating spectral differentiation matrices.
Otherwise, the default is to give the coefficients of the true definite
integral polynomial
  • Loading branch information
askhamwhat committed May 3, 2024
1 parent 4736fd2 commit 8012c8b
Show file tree
Hide file tree
Showing 4 changed files with 101 additions and 7 deletions.
23 changes: 23 additions & 0 deletions chunkie/+lege/dermat.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,23 @@
function mat = dermat(k,u,v)
%LEGE.DERMAT returns the spectral differentiation matrix on Legendre nodes
% of order k
%
% input:
% k - integer, number of Legendre nodes
%
% optional inputs:
% u - k x k matrix mapping values at legendre nodes to legendre series
% coefficients
% v - k x k matrix mapping legendre series coefficients to values at
% legendre nodes
%
% output:
% mat - spectral differentiation matrix on Legendre nodes
%
%see also LEGE.EXPS, LEGE.DERPOL

if nargin < 3
[~,~,u,v] = lege.exps(k);
end

mat = v(:,1:end-1)*lege.derpol(u);
19 changes: 14 additions & 5 deletions chunkie/+lege/intmat.m
Original file line number Diff line number Diff line change
@@ -1,15 +1,24 @@
function [aint,x,w] = intmat(n)
function [aint,u,v] = intmat(n,u,v)
%INTMAT returns the spectral integration matrix on n Gaussian nodes
% a transcription of part of the Rokhlin routine legeinmt
%
% input:
% n - the number of Gaussian nodes
%
% optional inputs:
% u - k x k matrix mapping values at legendre nodes to legendre series
% coefficients
% v - k x k matrix mapping legendre series coefficients to values at
% legendre nodes
%
% output:
%

[x,w,u,v] = lege.exps(n);
if nargin < 3
[~,~,u,v] = lege.exps(n);
end

coeffs = eye(n);
polints = lege.intpol(coeffs);
aint = v*polints(1:end-1,:)*u;
tmp = lege.intpol(u,'original');
aint = v*tmp(1:end-1,:);

end
34 changes: 32 additions & 2 deletions chunkie/+lege/intpol.m
Original file line number Diff line number Diff line change
@@ -1,12 +1,42 @@
function intcoeffs = intpol(coeffs)
function intcoeffs = intpol(coeffs,const_option)
%LEGE.INTPOL compute the coefficients of the indefinite
% integral of the polynomial with the coefficients given on
% input
%
% input:
% coeffs - array of Legendre series coefficients. if matrix, each column
% contains the Legendres series coefficients of a separate function
%
% optional input:
% const_option - string. default ('true'). If 'true', then the constant on
% output is such that the polynomial(s) described by intcoeffs on output
% take the value zero at -1. If 'original', then the constant on output
% is such that the polynomial(s) described by intcoeffs(1:end-1,:) on
% output take the value zero at -1. This is used for a spectral
% differentiation matrix taking point values to point values on a grid
% of the same order.
%
% output:
% intcoeffs - array of Legendre series coefficients for a definte integral
% of the input Legendre series
%

if nargin < 2
const_option = 'true';
end

sz = size(coeffs);
szint = sz; szint(1) = szint(1)+1;
intcoeffs = zeros(szint);

if strcmpi(const_option,'true')
ncc = szint(1);
elseif strcmpi(const_option,'original')
ncc = szint(1)-1;
else
error('LEGE.INTPOL: unknown option for constant');
end

for i = 2:sz(1)
j = i-1;
intcoeffs(i+1,:) = coeffs(i,:)/(2*j+1);
Expand All @@ -17,7 +47,7 @@

sss=-1;
dd = zeros(size(intcoeffs(end,:)));
for k = 2:(szint(1)-1)
for k = 2:ncc
dd=dd+intcoeffs(k,:)*sss;
sss=-sss;
end
Expand Down
32 changes: 32 additions & 0 deletions devtools/test/legeexpsunitTest.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
%

addpaths_loc();

k = 19;

[x,w,u,v] = lege.exps(k);

dmat = lege.dermat(k,u,v);
imat = lege.intmat(k,u,v);

pv = sin(x);
dpv_true = cos(x);
ipv_true = -cos(x)+cos(-1);

dpv = dmat*pv;
ipv = imat*pv;

assert(norm(dpv-dpv_true) < 1e-12)
assert(norm(ipv-ipv_true) < 1e-14)

cfs = randn(k,1);
cfsint1 = lege.intpol(cfs);
cfsintold = lege.intpol(cfs,'original');
fun = @(t) lege.exev(t,cfs);

for j = 1:k
itrue = integral(fun,-1,x(j));
icoefs = lege.exev(x(j),cfsint1);
assert(abs(itrue-icoefs)< 1e-14);
end

0 comments on commit 8012c8b

Please sign in to comment.