<br><br><div class="gmail_quote">On Mon, Dec 20, 2010 at 9:46 PM, Charles R Harris <span dir="ltr">&lt;<a href="mailto:charlesr.harris@gmail.com">charlesr.harris@gmail.com</a>&gt;</span> wrote:<br><blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
<br><br><div class="gmail_quote"><div><div></div><div class="h5">On Mon, Dec 20, 2010 at 9:17 PM, Charles R Harris <span dir="ltr">&lt;<a href="mailto:charlesr.harris@gmail.com" target="_blank">charlesr.harris@gmail.com</a>&gt;</span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
<br><br><div class="gmail_quote"><div><div></div><div>On Mon, Dec 20, 2010 at 8:55 PM, Eric Carlson <span dir="ltr">&lt;<a href="mailto:ecarlson@eng.ua.edu" target="_blank">ecarlson@eng.ua.edu</a>&gt;</span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
On 12/18/2010 12:35 PM, Pramod wrote:<br>
&gt; Matlab imple mentation :<br>
&gt; for N = 1:Nmax;<br>
&gt; [D,x] = cheb(N);<br>
&gt;<br>
&gt; How to impliment above (written in matlab ) chebshev  polynomial in<br>
&gt; python<br>
<br>
cheb is not a standard matlab function, but if this is it:<br>
<br>
function [D,x]=cheb(N)<br>
<br>
if N==0, D=0; x=1; return; end<br>
x=cos(pi*(0:N)/N)&#39;;<br>
c=[2; ones(N-1,1); 2].*(-1).^(0:N)&#39;;<br>
X=repmat(x,1,N+1);<br>
dX=X-X&#39;;<br>
D=(c*(1./c)&#39;)./(dX+eye(N+1)); % off diagonals<br>
D=D-diag(sum(D&#39;)); % diagonals<br>
<br>
<br>
Then a python version could be given by:<br>
<br>
from numpy import cos,pi,linspace,array,matrix,ones,hstack,eye,diag,tile<br>
<br>
def cheb(N):<br>
<br>
     if N==0:<br>
         D=0.0<br>
         x=1.0<br>
     else:<br>
         x=cos(pi*linspace(0,N,N+1)/N)<br>
         c=matrix(hstack(([2.],ones(N-1),[2.]))*(-1)**linspace(0,N,N+1)).T<br>
         X=tile(x,[N+1,1])<br>
         dX=(X-X.T).T<br>
         D=array(c*(1/c).T)/(dX+eye(N+1)); # off diagonals<br>
         D=D-diag(sum(D,axis=1)); # diagonals<br>
<br>
     return D,x<br>
<br>
(D,x)=cheb(3)<br>
print D<br>
(D,x)=cheb(4)<br>
print D<br>
<br>
<br>
I almost literally translated the matlab code, so I do not know how<br>
efficient it all is, and have given zero thought to better ways to do<br>
it. You need to be very careful with matrix and array types. Unless you<br>
KNOW you want a matrix, you probably want an array. One of the biggest<br>
transitions from matlab to python is learning to not worry about whether<br>
something is a column or row array, except when you are doing matrix<br>
multiplication.<br>
<br></blockquote></div></div><div><br>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.<br>


<br></div></div></blockquote></div></div><div><br>If so, the following brute force method will work.<br><br><span style="font-family: courier new,monospace;">In [1]: import numpy.polynomial as poly</span><br style="font-family: courier new,monospace;">

<br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">In [2]: def Cheb(N):</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">   ...:     x = poly.chebpts2(N+1)</span><br style="font-family: courier new,monospace;">

<span style="font-family: courier new,monospace;">   ...:     D = array([poly.Chebyshev.fit(x,c,N).deriv()(x) for c in eye(N+1)])</span><br style="font-family: courier new,monospace;"><span style="font-family: courier new,monospace;">   ...:     return D, x</span><br>

<br></div></div></blockquote><div><br>Oops, D needs a transpose.<br><br>Chuck <br></div></div>