Code covered by the BSD License

Chebfun V4

Chebfun Team (view profile)

30 Apr 2009 (Updated )

Numerical computation with functions instead of numbers.

Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

jacpoly(N,a,b,d,flag)
function Q = jacpoly(N,a,b,d,flag)
%JACPOLY Jacobi polynomials.
%   P = JACPOLY(N,ALPHA,BETA) computes a chebfun of the Jacobi polynomial
%   of degree N with parameters ALPHA. BETA on the interval [-1,1].
%   N can be a vector of integers.
%
%   P = JACPOLY(N,ALPHA,BETA,D) computes the Jacobi polynomials as above,
%   but on the interval given by the domain D, which must be bounded.
%
%   JACPOLY(N,ALPHA,BETA,D,'norm') or JACPOLY(N,ALPHA,BETA,'norm') will
%   normalise so that P(1) or P(D(2)) = 1 ;
%

% Copyright 2011 by The University of Oxford and The Chebfun Developers.
% See http://www.maths.ox.ac.uk/chebfun/ for Chebfun information.

Nmax = max(N);
x = chebpts(Nmax+1);

if nargin < 3, error('CHEBFUN:jacpoly:inputs','jacpoly requires at least 3 inputs'); end
if nargin < 4, d = [-1,1]; end
if nargin < 5, flag = 0; end
if isa(d,'char'), flag = d; d = [-1,1]; end
if isa(d,'domain'), d = d.ends; end

apb = a+b;
P = zeros(Nmax+1);
P(:,1) = 1;
P(:,2) = 0.5*(2*(a+1) + (apb+2)*(x-1));
for k = 2:Nmax
k2 = 2*k;
k2apb = k2+apb;
q1 =  k2*(k + apb)*(k2apb - 2);
q2 = (k2apb - 1)*(a*a - b*b);
q3 = (k2apb - 2)*(k2apb - 1)*k2apb;
q4 =  2*(k + a - 1)*(k + b - 1)*k2apb;
P(:,k+1) = ((q2+q3*x).*P(:,k)-q4*P(:,k-1))/q1;
end

Q = chebfun;
for k = 1:length(N)
xk = chebpts(N(k)+1,d);
Pk = P(:,N(k)+1);
if flag
Pk = Pk/Pk(end);
end
Q(:,k) = chebfun(bary(xk,Pk),d);
end