[SciPy-User] Chebeshev Polynomial implimentation python
Eric Carlson
ecarlson at eng.ua.edu
Mon Dec 20 22:55:09 EST 2010
On 12/18/2010 12:35 PM, Pramod wrote:
> Matlab imple mentation :
> for N = 1:Nmax;
> [D,x] = cheb(N);
>
> How to impliment above (written in matlab ) chebshev polynomial in
> python
cheb is not a standard matlab function, but if this is it:
function [D,x]=cheb(N)
if N==0, D=0; x=1; return; end
x=cos(pi*(0:N)/N)';
c=[2; ones(N-1,1); 2].*(-1).^(0:N)';
X=repmat(x,1,N+1);
dX=X-X';
D=(c*(1./c)')./(dX+eye(N+1)); % off diagonals
D=D-diag(sum(D')); % diagonals
Then a python version could be given by:
from numpy import cos,pi,linspace,array,matrix,ones,hstack,eye,diag,tile
def cheb(N):
if N==0:
D=0.0
x=1.0
else:
x=cos(pi*linspace(0,N,N+1)/N)
c=matrix(hstack(([2.],ones(N-1),[2.]))*(-1)**linspace(0,N,N+1)).T
X=tile(x,[N+1,1])
dX=(X-X.T).T
D=array(c*(1/c).T)/(dX+eye(N+1)); # off diagonals
D=D-diag(sum(D,axis=1)); # diagonals
return D,x
(D,x)=cheb(3)
print D
(D,x)=cheb(4)
print D
I almost literally translated the matlab code, so I do not know how
efficient it all is, and have given zero thought to better ways to do
it. You need to be very careful with matrix and array types. Unless you
KNOW you want a matrix, you probably want an array. One of the biggest
transitions from matlab to python is learning to not worry about whether
something is a column or row array, except when you are doing matrix
multiplication.
HTH
More information about the SciPy-User
mailing list