High-degree, zeroth-order Legendre Polynomials (computation speed)

127 views
Skip to first unread message

CW Edwards

unread,
Apr 29, 2021, 12:07:31 PM4/29/21
to chebfun-users
I stumbled upon the Chebfun software package and supporting scientific papers while in search for a fast method of computing (i.e. approximating) high-degree, zeroth-order Legendre polynomials in MATLAB. I am currently using MATLAB's built-in Legendre function (legendre.m), which is proving to be a bottleneck in my code.

I believe this is due to the fact that, by default, MATLAB's Legendre function computes and returns all integer-valued orders of Legendre functions up to the specified degree. That is,

P = legendre(N,X) computes the associated Legendre functions 
    of degree N and order M = 0, 1, ..., N, evaluated for each element
    of X.  N must be a scalar integer and X must contain real values
    between -1 <= X <= 1.  
My hope was that by using the Chebfun package to compute only the zeroth-order polynomials, I would gain some efficiency in the code. As a preliminary test, I simply timed the two functions (i.e. MATLAB's legendre.m versus CHEBFUN's legpoly.m) while computing a Legendre polynomial of degree 100 over a densely sampled domain of [-1,1]. Surprisingly, the CHEBFUN function was slower than MATLAB's function by an order of magnitude even though the latter is computing all of the lower order functions (i.e. M = 0, 1, …, 100). 

I assume the results I've obtained here are due to my ignorance of how CHEBFUN works. Perhaps legpoly.m is not the actual function that I should be using? 

For clarity, I'd ultimately like to approximate the value of these high-degree, zeroth-order Legendre polynomials for a large matrix of input arguments. 

Any assistance that can be provided here is greatly appreciated. 

Jacob Rus

unread,
Apr 29, 2021, 12:30:42 PM4/29/21
to CW Edwards, chebfun-users
Hi CW, I don't see the same kind of timing you report.

>> addpath('path/to/chebfun')
>> X = rand(1e6);
>> N = 100;
>> tic; legendre(N, X); toc
Elapsed time is 7.6 seconds
>> tic; L = legpoly(N); L(X); toc;
Elapsed time is 0.22 seconds

Cheers,
Jacob

Jacob Rus

unread,
Apr 29, 2021, 12:33:07 PM4/29/21
to CW Edwards, chebfun-users
>> X = rand(1e6);

erm, typo alert. Should be `X = rand(1e6, 1);`

–Jacob

Austin, Anthony (CIV)

unread,
Apr 29, 2021, 12:44:16 PM4/29/21
to CW Edwards, chebfun-users
Hi CW,

Performance can be a tricky thing to troubleshoot since it's very
machine- and code-dependent. Would it be possible for you to share the
code that does your performance test?

For what it's worth, I would evaluate the degree-100 Legendre polynomial
using Chebfun like this:

x = ... % [Vector of points at which to evaluate]
P = legpoly(100);
y = P(x);

Here's a short script that compares this approach to MATLAB's built-in
legendre():

-----[BEGIN CODE]-----

n = 100; % Polynomial degree
N = 10000; % No. of evaluation points
Ntrials = 5;

x = linspace(-1, 1, N).';

y1 = legendre(n, x);

tic();
for (i = 1:1:Ntrials)
y1 = legendre(n, x);
end
t1 = toc()/Ntrials;
fprintf('MATLAB (legendre): % .2e\n', t1);

Pn = legpoly(n);
y2 = Pn(x);

tic();
for (i = 1:1:Ntrials)
Pn = legpoly(n);
y2 = Pn(x);
end
t2 = toc()/Ntrials;
fprintf('Chebfun (build + eval): % .2e\n', t2);

Pn = legpoly(n);
y3 = Pn(x);

tic();
for (i = 1:1:Ntrials)
y3 = Pn(x);
end
t3 = toc()/Ntrials;
fprintf('Chebfun (eval only): % .2e\n', t3);

-----[END CODE]-----

A typical run on my machine prints the following:

MATLAB (legendre): 4.11e-02
Chebfun (build + eval): 7.28e-03
Chebfun (eval only): 4.22e-03

That is, contrary to your results, I find Chebfun to be 10x faster than
legendre() for this evaluation. (If you build the chebfun before every
evaluation, it's only 5x faster, but you can/should avoid this in
practice.)


Regards,

Anthony P. Austin
anthony...@nps.edu
> --
> You received this message because you are subscribed to the Google Groups
> "chebfun-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email
> to chebfun-user...@googlegroups.com.
> To view this discussion on the web, visithttps://groups.google.com/d/msgid/chebfun-users/b800320a-5959-41be-8acb-576b4972
> 5becn%40googlegroups.com.

CW Edwards

unread,
Apr 29, 2021, 1:07:28 PM4/29/21
to chebfun-users
Jacob and Anthony,

Thank you for the replies! I've included the code that I am currently using to time the two methods below. 

The timing results that are being reported by MATLAB are:

Elapsed time is 0.002763 seconds. (MATLAB's legendre.m)
Elapsed time is 0.070896 seconds. (CHEBFUN's legpoly.m)

-----[BEGIN CODE]-----

X = -1:.01:1;
N = [0 1 2 3 4 5 6 100];

% Compute Legendre Polynomials using MATLAB's legendre.m 
P_N = zeros(length(N),length(X));
tic
for n = 1:length(N)
    P = legendre(N(n),X);
    if n > 0
        P = squeeze(P(1,:,:));
    end
    P_N(n,:) = P;
end
toc

% Compute Legendre Polynomials using CHEBFUN's legpoly.m
tic
P_N_CHEBFUN = legpoly(N,X);
toc

% Plot output from MATLAB's legendre.m 
fig_h1 = figure();
plot(X,P_N);
legend('$P_0$','$P_1$','$P_2$','$P_3$','$P_4$','$P_5$','$P_6$','$P_{100}$','location','southeast','interpreter','latex');
title('Legendre Polynomials (MATLAB)');

% Plot output from CHEBFUN's legpoly.m 
fig_h2 = figure();
plot(P_N_CHEBFUN); 
legend('$P_0$','$P_1$','$P_2$','$P_3$','$P_4$','$P_5$','$P_6$','$P_{100}$','location','southeast','interpreter','latex');
title('Legendre Polynomials (CHEBFUN)');

-----[END CODE]-----

CW

Austin, Anthony (CIV)

unread,
Apr 29, 2021, 1:15:54 PM4/29/21
to CW Edwards, chebfun-users
> P_N_CHEBFUN = legpoly(N,X);

This line does not do what you want. It just builds a piecewise chebfun
for the degree-N Legendre polynomial over the domain specified in X
(interpreting the entries in X as breakpoints).

The right way to do your evaluation in Chebfun is

P = legpoly(100);
P_N_CHEBFUN = P(X);

As my code showed, that ought to be faster than MATLAB's legendre(), at
least if the number of evaluation points is large enough.


Regards,

Anthony P. Austin
anthony...@nps.edu

On Thu, 29 Apr 2021, CW Edwards wrote:

> --
> You received this message because you are subscribed to the Google Groups
> "chebfun-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email
> to chebfun-user...@googlegroups.com.
> To view this discussion on the web, visithttps://groups.google.com/d/msgid/chebfun-users/cf5bc96d-0f6d-4d10-a6c0-2f3dcb46
> e07cn%40googlegroups.com.
>
>

CW Edwards

unread,
Apr 29, 2021, 1:28:48 PM4/29/21
to chebfun-users
Anthony,

Thanks again! 

Is it possible to use legpoly.m as you've suggested if N remains a vector?

If I set N = 100 (i.e. a scalar) in my previous code and use legpoly as suggested, I do obtain a more favorable result with both methods resulting in a similar execution time (though legendre.m is still reported as being slightly faster).

Conversely, if I let N remain a vector as in the original code and execute the following:

P_N_CHEBFUN = legpoly(N);
Y = P_N_CHEBFUN(X);

The resulting data structure Y is a row vector of size  1x1608, which I'm not certain how to interpret. 

To obtain my desired result would I need to use legpoly in a loop similar to how I've implemented use of legendre.m in my original code?

All the best,

CW

CW Edwards

unread,
Apr 29, 2021, 1:47:49 PM4/29/21
to chebfun-users
To better illustrate my last question I've included updated code below.

Elapsed time is 0.002220 seconds. (MATLAB's legendre.m)
Elapsed time is 0.017406 seconds. (CHEBFUN's legpoly.m)

As can be seen, I am still getting a shorter execution time (by a order of magnitude) with legendre.m, if I'm required to execute both functions in a loop to get results for Legendre polynomials of multiple degrees.

The major concern here being that if legendre.m is faster for a relatively small set of inputs, as specified in the vector X below, the gap may widen even more for a large matrix of inputs. 

As before, any clarity that can be added here is greatly appreciated. 

-----[BEGIN CODE]-----

X = -1:.01:1;
N = [0 1 2 3 4 5 6 100];

% Compute Legendre Polynomials using MATLAB's legendre.m 
P_N = zeros(length(N),length(X));
tic
for n = 1:length(N)
    P = legendre(N(n),X);
    if n > 0
        P = squeeze(P(1,:,:));
    end
    P_N(n,:) = P;
end
toc

% Compute Legendre Polynomials using CHEBFUN's legpoly.m
P_N_CHEBFUN = zeros(length(N),length(X));
tic
for n = 1:length(N)
    P = legpoly(N(n));
    P_N_CHEBFUN(n,:) = P(X);
end
toc

% Plot output from MATLAB's legendre.m 
fig_h1 = figure();
plot(X,P_N);
legend('$P_0$','$P_1$','$P_2$','$P_3$','$P_4$','$P_5$','$P_6$','$P_{100}$','location','southeast','interpreter','latex');
title('Legendre Polynomials (MATLAB)');

% Plot output from CHEBFUN's legpoly.m 
fig_h2 = figure();
plot(X,P_N_CHEBFUN); 
legend('$P_0$','$P_1$','$P_2$','$P_3$','$P_4$','$P_5$','$P_6$','$P_{100}$','location','southeast','interpreter','latex');
title('Legendre Polynomials (CHEBFUN)');

-----[END CODE]-----

Austin, Anthony (CIV)

unread,
Apr 29, 2021, 2:12:06 PM4/29/21
to CW Edwards, chebfun-users
> Is it possible to use legpoly.m as you've suggested if N remains a
> vector?

If N is a vector, legpoly(N) builds a chebfun quasimatrix whose columns
are the Legendre polynomials corresponding to the entries in N. For
instance, after

N = [0 1 2 100];
P = legpoly(N);

P(:, 1) will be the Legendre polynomial of degree N(1) = 0, P(:, 2) will
be the Legendre polynomial of degree N(2) = 1, etc.

Point evaluations with this construct work as follows. If x is a column
vector, then after

Y = P(x);

The kth column Y(:, k) of Y will contain the values of P(:, k) at x.
If you only want to evaluate one polynomial, say, the one whose degree
is N(3), you can do

y = P(x, 3);

In general, I would avoid evaluating quasimatrices at row vectors
because the results, while sensible, aren't usually convenient to work
with. It's almost always easier to convert the input to a column vector
and then reshape the output as needed.

> To obtain my desired result would I need to use legpoly in a loop
> similar to how I've implemented use of legendre.m in my original code?

Sometimes doing the vectorized evaluation like the above is faster.
Somteimes, it's faster to loop through the desired degree and build
individual chebfuns for each polynomial. You have to experiment to see.

> If I set N = 100 (i.e. a scalar) in my previous code and use legpoly
> as suggested, I do obtain a more favorable result with both methods
> resulting in a similar execution time (though legendre.m is still
> reported as being slightly faster).

This probably happens because your evaluation grid is pretty small: 201
points may not be enough for Chebfun to make a difference in your
environment.

My own code with the no. of evaluation points set to 201 prints

MATLAB (legendre): 1.26e-03
Chebfun (build + eval): 3.75e-03
Chebfun (eval only): 2.42e-04

Chebfun is still about 5x faster in evaluation; it's only if you include
the time to build the chebfun that it becomes slower. If your
application involves doing a large number of evaluations, you should
build the chebfun just once and reuse it each time so you can amortize
the cost of the construction over all the evaluations.

Austin, Anthony (CIV)

unread,
Apr 29, 2021, 2:12:13 PM4/29/21
to CW Edwards, chebfun-users
> The major concern here being that if legendre.m is faster for a relatively small
> set of inputs, as specified in the vector X below, the gap may widen even more
> for a large matrix of inputs. 

No! It is *exactly* the opposite! Chebfun (properly used) will be
*faster* than legendre() for large inputs, while legendre() will
typically be faster for small inputs, since for small inputs the
overhead of chebfun construction is too large for the faster evaluation
to be worth it.

Try this:

-----[BEGIN CODE]-----

% A few different grid sizes---try them!
%X = -1:.01:1;
%X = -1:.001:1;
%X = -1:.0001:1;
X = -1:.00001:1;

N = [0 1 2 3 4 5 6 100];

% Compute Legendre Polynomials using MATLAB's legendre.m
P_N = zeros(length(N),length(X));
tic
for n = 1:length(N)
P = legendre(N(n),X);
if n > 0
P = squeeze(P(1,:,:));
end
P_N(n,:) = P;
end
toc

% Compute Legendre polynomials with Chebfun's legpoly(), evaluating in a
% vectorized fashion.
tic
P = legpoly(N);
P_N_CHEBFUN = P(X.').';
toc

% Compute Legendre polynomials with Chebfun's legpoly(), building
% individual chebfuns for each degree.
P_N_CHEBFUN_2 = zeros(length(N), length(X));
tic
for n = 1:length(N)
P = legpoly(N(n));
P_N_CHEBFUN_2(n, :) = P(X).';
end
toc

% Plot output from MATLAB's legendre.m
fig_h1 = figure(1);
plot(X,P_N);
legend('$P_0$','$P_1$','$P_2$','$P_3$','$P_4$','$P_5$','$P_6$','$P_{100}$', 'location','southeast','interpreter','latex');
title('Legendre Polynomials (MATLAB)');

% Plot output from CHEBFUN's legpoly.m
fig_h2 = figure(2);
plot(P_N_CHEBFUN.');
legend('$P_0$','$P_1$','$P_2$','$P_3$','$P_4$','$P_5$','$P_6$','$P_{100}$','location','southeast','interpreter','latex');
title('Legendre Polynomials (CHEBFUN)');

% Plot output from CHEBFUN's legpoly.m
fig_h3 = figure(3);
plot(P_N_CHEBFUN_2.');
legend('$P_0$','$P_1$','$P_2$','$P_3$','$P_4$','$P_5$','$P_6$','$P_{100}$','location','southeast','interpreter','latex');
title('Legendre Polynomials (CHEBFUN)');

-----[END CODE]-----

With X = -1:0.01:1 selected, I get

Elapsed time is 0.003227 seconds.
Elapsed time is 0.008599 seconds.
Elapsed time is 0.017367 seconds.

So MATLAB's legendre() is faster here. With X = -1:0.00001:1 selected, I get

Elapsed time is 0.906905 seconds.
Elapsed time is 0.402918 seconds.
Elapsed time is 0.077085 seconds.

So Chebfun is faster here, and it is better to build individual chebfuns
than to do a vectorized evaluation with a quasimatrix.

Nick Hale

unread,
Apr 29, 2021, 3:07:28 PM4/29/21
to chebfun-users
Hi all

I'm a little late to the party. My inputs:

 * I don't know why the vectorised version of Clenshaw is so slow here. 

 * If speed is the goal, and you don't care about syntax (which if you're
   considering the legendre->squeeze option, then you don't), then you
   will be hard pressed to beat this (which avoids constructor overhead):

tic
I = eye(max(N+1));
c_leg = I(:,N+1);
c_cheb = leg2cheb(c);
Y = zeros(length(X), size(c,2));
for k = 1:size(c,2)
    Y(:,k) = chebtech.clenshaw(X',c_cheb(:,k));
end
toc

   (The for loop could be avoided and 
      Y = chebtech.clenshaw(X',c_cheb); 
   would suffice if vectorised Clenshaw wasn't so slow!)

N

PS. I don't consider 100 to be particularly high degree! Try N = 1000.

Nick Hale

unread,
Apr 29, 2021, 3:12:35 PM4/29/21
to chebfun-users
Google Groups is archaic and I can't edit my post. Correction to my code here:

tic
I = eye(max(N+1));
c_leg = I(:,N+1);
c_cheb = leg2cheb(c_leg);
Y = zeros(length(X), size(c_cheb,2));
for k = 1:size(c_cheb,2)
    Y(:,k) = chebtech.clenshaw(X',c_cheb(:,k));
end
toc

Apologies for spam.

N

Austin, Anthony (CIV)

unread,
Apr 29, 2021, 3:24:12 PM4/29/21
to Nick Hale, chebfun-users
>  * I don't know why the vectorised version of Clenshaw is so slow here. 

I'm pretty sure the reason the quasimatrix eval. / vectorized Clenshaw
is slow here vs. evaluating each chebfun individually is due mainly to
the discrepancy in the degrees. When we build legpoly([0 1 2 3 4 5 6
100]), all 8 columns have to be represented with length 101, which is a
bit of a waste as far as the first 7 columns are concerned. If I force
an array-of-chebfuns representation with cheb2quasi, the time comes way
down:

-----[BEGIN CODE]-----
x = (-1:.00001:1).';
N = [0 1 2 3 4 5 6 100];

P = legpoly(N);
Y = P(x);

tic();
Y = P(x);
toc();

P = cheb2quasi(legpoly(N));
Y = P(x);

tic();
Y = P(x);
toc();

P = cell(size(N));
for (n = 1:1:length(N))
P{n} = legpoly(N(n));
end

tic();
for n = 1:length(N)
Y = P{n}(x);
end
toc();
-----[END CODE]-----

A typical run of this on my machine prints

Elapsed time is 0.598007 seconds.
Elapsed time is 0.084444 seconds.
Elapsed time is 0.028718 seconds.

I guess the remaining factor of 3 is due to whatever overhead there is
in the array-of-chebfuns representation versus individual ones.

Austin, Anthony (CIV)

unread,
Apr 29, 2021, 3:55:57 PM4/29/21
to Nick Hale, chebfun-users
Greetings! It's been a while. How have you been, waiting out the
pandemic in Stellenbosch?

> tic
> I = eye(max(N+1));
> c_leg = I(:,N+1);
> c_cheb = leg2cheb(c_leg);
> Y = zeros(length(X), size(c_cheb,2));
> for k = 1:size(c_cheb,2)
>     Y(:,k) = chebtech.clenshaw(X',c_cheb(:,k));
> end
> toc

I actually measure this to still be slower than building the chebfuns
individually:

-----[BEGIN CODE]-----
x = (-1:.00001:1).';
N = [0 1 2 3 4 5 6 100];

% 1. Vectorized eval.
P = legpoly(N);
Y = P(x);

tic();
Y = P(x);
toc();

% 2. Individual chebfuns.
P = cell(size(N));
for (n = 1:1:length(N))
P{n} = legpoly(N(n));
end

Y = zeros(length(x), length(N));

tic();
for (n = 1:1:length(N))
Y(:, n) = P{n}(x);
end
toc();

% 3. Plumbing functions.
I = eye(max(N + 1));
c_leg = I(:, N + 1);
c_cheb = leg2cheb(c_leg);

Y = zeros(length(x), length(N));

tic();
for (n = 1:1:length(N))
Y(:, n) = chebtech.clenshaw(x, c_cheb(:, n));
end
toc();
-----[END CODE]-----

For me, this prints:

Elapsed time is 0.485624 seconds.
Elapsed time is 0.037628 seconds.
Elapsed time is 0.077818 seconds.

The reason is that the vectorized leg2cheb() call still makes some
unnecessarily long vectors for the low degrees here. If I add

-----[BEGIN CODE]-----
% 4. Plumbing functions, building coeffs vectors one a time:
c_cheb = cell(size(N));
for (n = 1:1:length(N))
I = eye(max(N(n) + 1));
c_leg = I(:, N(n) + 1);
c_cheb{n} = leg2cheb(c_leg);
end

Y = zeros(length(x), length(N));

tic();
for (n = 1:1:length(N))
Y(:, n) = chebtech.clenshaw(x, c_cheb{n});
end
toc();
-----[END CODE]-----

and run again, now I get a fourth line:

Elapsed time is 0.014742 seconds.

Much better.

Actually, your original code shows that the performance problem with
vectorized Clenshaw is subtler than I thought. Whatever it is, it's not
just that some of the columns were excessively long. (Your original
code shows that the penalty for that should be only about a factor of 2
slower than building individual chebfuns, while the quasimatrix eval. is
slower by more like a factor of 10...)


Regards,

Anthony P. Austin
anthony...@nps.edu

> --
> You received this message because you are subscribed to the Google Groups
> "chebfun-users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email
> to chebfun-user...@googlegroups.com.
> To view this discussion on the web, visithttps://groups.google.com/d/msgid/chebfun-users/a553f3b1-c9d3-4c4b-9e0f-e3830274
> da7cn%40googlegroups.com.
>
>

CW Edwards

unread,
Apr 29, 2021, 7:38:33 PM4/29/21
to chebfun-users
All,

Catching up here. Let me first say thank you to everyone for the very informative responses. 

The N vector used to specify the degrees of the Legendre Polynomials in my original code was simply intended to provide a range of different polynomials that could be plotted for visual verification.

In light of Anthony's comments regarding the sparseness of this vector, I have rewritten my example code to be more in-line with my actual application. 

In particular, the updated code specifies N = 0, 1, ..., 100 and evaluates the approximations of the Legendre Polynomials for a 200 x 200 matrix of input arguments. 

Following suit to Anthony's last post, I perform the calculation using three different methods:

(1.) Computing individual CHEBFUNS and calling legpoly.m
(2.) Computing individual CHEBFUNS and calling the Clenshaw method directly (thanks to Nick for this contribution)
(3.) Vectorizing CHEBFUNS and calling the Clenshaw method directly

The execution times that MATLAB reports for these three methods are:

Elapsed time is 0.138154 seconds. (Method 1)
Elapsed time is 0.096231 seconds. (Method 2)
Elapsed time is 0.147856 seconds. (Method 3)

Taking this into account, it would seem that Method 2 is the fastest, though only slightly faster than Method 1.

I do find it odd that the vectorized version (i.e. Method 3) performs worse given that we are executing this code in MATLAB. 

If anyone has thoughts on how to speed up the calculation even further I would be grateful to hear them.

Here's the updated example code:

-----[BEGIN CODE]-----

numPTerms = 100; % specify highest degree of zeroth-order Legendre functions to approximate
numInputs = 200; % specify number of terms in each dimension of input matrix
Lparg = rand([numInputs,numInputs]); % create input matrix (i.e. argument to Legendre Polynomials)

%% Compute Legendre Polynomials by computing individual CHEBFUNS and calling legpoly.m
P = cell(numPTerms+1);
for nn = 1:numPTerms+1
    P{nn} = legpoly(nn-1);
end

% Compute Legendre polynomials from CHEBFUN cell array
Pn_MATLAB = zeros(numInputs,numInputs,numPTerms+1); % preallocate memory for 3D matrix
tic
for nn = 1:numPTerms+1
   Pn_MATLAB(:,:,nn) = P{nn}(Lparg);
end
toc

%% Compute Legendre Polynomials by computing individual CHEBFUNS and calling Clenshaw method directly
N = 0:1:numPTerms;
c_cheb = cell(size(N));
for n = 1:1:length(N)
    I = eye(max(N(n) + 1));
    c_leg = I(:, N(n) + 1);
    c_cheb{n} = leg2cheb(c_leg);
end

Pn_CHEBFUN = zeros(numInputs,numInputs,numPTerms+1); % preallocate memory for 3D matrix
Lparg_vec = Lparg(:);
tic
for nn = 1:numPTerms+1
     Pn_CHEBFUN(:,:,nn) = reshape(chebtech.clenshaw(Lparg_vec, c_cheb{nn}),numInputs,numInputs);
end
toc

%% Compute Legendre Polynomials by vectorizing CHEBFUNS and calling Clenshaw method directly
N = 0:1:numPTerms;
I = eye(max(N + 1));
c_leg = I(:, N + 1);
c_cheb_2 = leg2cheb(c_leg);

Pn_CHEBFUN_2 = zeros(numInputs,numInputs,numPTerms+1); % preallocate memory for 3D matrix
Lparg_vec = Lparg(:);
tic
for nn = 1:numPTerms+1
     Pn_CHEBFUN_2(:,:,nn) = reshape(chebtech.clenshaw(Lparg_vec, c_cheb_2(:, nn)),numInputs,numInputs);
end
toc

-----[END CODE]-----

Jacob Rus

unread,
Apr 29, 2021, 9:11:54 PM4/29/21
to CW Edwards, chebfun-users
Hi CW,

If you really need maximum performance of Clenshaw evaluation, a
simple SIMD implementation in C is 1–2 orders of magnitude faster than
the Matlab implementation.

Here is some code:
https://github.com/alkashani/numcrunch/blob/master/src/clenshaw_vector.c

cf. https://github.com/chebfun/chebfun/issues/1582

Cheers,
Jacob

Nick Hale

unread,
Apr 30, 2021, 3:18:43 AM4/30/21
to chebfun-users
If you need all polynomials of degree less than or equal to N, then you can simply use the recurrence relation.

Another approach is to use the weighted QR implemented in legpoly (which is better for maintaining orthogonality).

NH

Ian Ajzenszmidt

unread,
May 1, 2021, 10:51:48 AM5/1/21
to CW Edwards, chebfun-users
Result of a scholar.google.com search on 

Ian
Scholar

Generalized energy condensation theory

F Rahnema, S Douglass, B Forget - Nuclear Science and …, 2008 - Taylor & Francis
… and there- fore allows for a high degree of flexibility in application … of the general theory. The additional moment equations coupled to the zeroth order represent correction terms to the flat flux ~in energy … Pl~m ! l'th order Legendre polynomials …
Cited by 37 Related articles

Experimental analysis of shape deformation of evaporating droplet using Legendre polynomials

A Sanyal, S Basu, R Kumar - Physics Letters A, 2014 - Elsevier
… k z , k r – axial and radial wave numbers, J 0 – zeroth order Bessel function … droplet shape read from the images and the reconstructed shape using Legendre polynomial expansion (Fig … feature into our analysis except in case of kerosene which shows high degree of deformation …
Cited by 12 Related articles

FFT-based high-performance spherical harmonic transformation

C Gruber, P Novák, J Sebera - Studia Geophysica et Geodaetica, 2011 - Springer
… coefficients with the Fourier constants of the associated Legendre functions, and (3) a subsequent 2 … The simplification of Eq.(9) is feasible for the real-valued function V. In the … spherical harmonic coefficients is an over-determined system, except for the zeroth order (see Sneeuw …
Cited by 32 Related articles

Computer-Aided Calibration of a Fine Autorod

CE Cohn, RA Karam, JE Marshall… - Nuclear Applications …, 1969 - Taylor & Francis
… is important that its cali- bration be known to a high degree of accuracy … Its inclusion makes redundant the zeroth-order Legendre and Gram polynomials, which are also constants … diagonal if the conditions for orthogonality are perfectly met.) Legendre polynomials are orthogonal …
Cited by 5 Related articles

Computational aspects of a code to study rotating turbulent convection in spherical shells

TC Clune, JR Elliott, MS Miesch, J Toomre… - Parallel Computing, 1999 - Elsevier
… parts of the perturbations ρ and P, and T is obtained from the zeroth-order equation of … may easily be accomplished in spectral space by making use of the appropriate Legendre and Fourier … S l m (r,t)) in terms of Chebyshev polynomials (chosen for their high degree of accuracy …
Cited by 209 Related articles

Modeling real shim fields for very high degree (and order) B0 shimming of the human brain at 9.4 T

P Chang, S Nassirpour… - Magnetic resonance in …, 2018 - Wiley Online Library
… Spherical harmonic shim systems are based on analytical models (from the Legendre polynomials), and vendor … by each of the shim coils are identical to the desired spherical harmonic function 15, 16 … This has been done previously for very highdegree (or very high‐order) shim …
Cited by 24 Related articles

Computing many-body wave functions with guaranteed precision: The first-order Møller-Plesset wave function for the ground state of helium atom

FA Bischoff, RJ Harrison, EF Valeev - The Journal of chemical …, 2012 - aip.scitation.org
… element representation with the number of dimensions, a brute-force computation of two … are typically used to eliminate the high-frequency components of the wave function … Wavelets for density-functional theory and post-density-functional theory calculations,” in Theoretical and …
Cited by 44 Related articles

Massively parallel method of characteristics neutron transport calculation with anisotropic scattering treatment on GPUs

N Choi, J Kang, H Joo - … on High Performance Computing in Asia-Pacific …, 2018 - dl.acm.org
… Nuclear power reactors are characterized by high degree of radial heterogeneity and small axial heterogenity due to … Pl is the ordinary Legendre polynomial of order l, α is the azimuthal angle, and μ … is the l-th order Legendre angular expanded scattering cross section, and μs is …
Cited by 10 Related articles

Efficient interpolation of high-fidelity geopotentials

N Arora, RP Russell - Journal of Guidance, Control, and Dynamics, 2016 - arc.aiaa.org
… The cubed-sphere model is discontinuous even to the zeroth order across the shell boundaries, although the … latitude of the position vector, respectively; and P 2 is the second-degree Legendre polynomial … In the current study, simple polynomials are chosen to minimize runtime …
Cited by 13 Related articles

Collimated light sources in the diffusion approximation

T Spott, LO Svaasand - Applied optics, 2000 - osapublishing.org
… It should describe the source as having a high degree of isotropy and at the same time should mimic the effects of … A.2 shows that the Legendre expansion coefficients are given by qk z 2 1 z . They are thus equal at … The zeroth-order moment of the photon distribution is given by …
Cited by 94 Related articles

--
You received this message because you are subscribed to the Google Groups "chebfun-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to chebfun-user...@googlegroups.com.
Reply all
Reply to author
Forward
0 new messages