[SciPy-User] Chebeshev Polynomial implimentation python

[Charles R Harris]

On Mon, Dec 20, 2010 at 9:46 PM, Charles R Harris <charlesr.harris at gmail.com
> wrote:

>
>
> On Mon, Dec 20, 2010 at 9:17 PM, Charles R Harris <
> charlesr.harris at gmail.com> wrote:
>
>>
>>
>> On Mon, Dec 20, 2010 at 8:55 PM, Eric Carlson <ecarlson at 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.
>>
>>
> If so, the following brute force method will work.
>
> In [1]: import numpy.polynomial as poly
>
> In [2]: def Cheb(N):
>    ...:     x = poly.chebpts2(N+1)
>    ...:     D = array([poly.Chebyshev.fit(x,c,N).deriv()(x) for c in
> eye(N+1)])
>    ...:     return D, x
>
>
Oops, D needs a transpose.

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