[SciPy-User] Chebeshev Polynomial implimentation python

Charles R Harris charlesr.harris@gmail....
Mon Dec 20 22:17:00 CST 2010


On Mon, Dec 20, 2010 at 8:55 PM, Eric Carlson <ecarlson@eng.ua.edu> wrote:

> 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.
>
>
As a guess, D is the differentiation operator for a function sampled at the
Chebyshev points of the second kind. So perhaps this is intended as a method
of solution for a boundary value problem.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20101220/b86bd1d5/attachment.html 


More information about the SciPy-User mailing list