[SciPy-User] Chebeshev Polynomial implimentation python

[Charles R Harris]

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

However, I have a module for this type of use case which I've attached.  I'm
not sure what shape it's in, it's old. The relevant function would be
modified_derivative

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20101220/a6693790/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: chebyshev.py
Type: text/x-python
Size: 18606 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20101220/a6693790/attachment.py>