home *** CD-ROM | disk | FTP | other *** search
/ Power Programming / powerprogramming1994.iso / progtool / microcrn / issue_40.arc / DAIMS.ARC / CHEB_VEC.CPP < prev    next >
Text File  |  1988-02-10  |  3KB  |  146 lines
  1. #include <math.h>
  2. #include "matrix.hxx"
  3. #include "vimatrix.hxx"
  4. #include "Cheb_vector.hxx"
  5.  
  6. double Cheb_vector::dx() { return  2/(size()-1);} 
  7.                           /* phys space range is from -1 to 1 */
  8. /*
  9. -*++ Cheb_vector::physical(): transforms a Chebyshev vector into physical space
  10. ** 
  11. ** (*++ history: 
  12. **     17 Dec 87    Bruce Eckel    Creation date
  13. ** ++*)
  14. ** 
  15. ** (*++ detailed: 
  16. ** ++*)
  17. */
  18.  
  19. phys_vector & Cheb_vector::physical()
  20. {
  21.     int mode = 1;
  22.     phys_vector & physical_rep = *new phys_vector(size(),-1);
  23.     double fact = M_PI * mode / size();
  24.     for (int i = 0; i < size(); i++) {
  25.     double X = cos( ((double)i + 0.5) * fact);
  26.     double summation =  0;
  27.     for (int n = 0; n < size(); n++)
  28.         summation += (*this)[n] * T(n,X);
  29.     physical_rep[i] = summation;
  30.     }
  31.     return physical_rep;
  32. }
  33.  
  34. /* 
  35. -*++ Cheb_vector::prime(): Returns the first derivative of a Cheb_vector
  36. ** 
  37. ** (*++ history: 
  38. **     17 Dec 87    Bruce Eckel    Creation date
  39. ** ++*)
  40. ** 
  41. ** (*++ detailed: 
  42. ** ++*)
  43. */
  44.  
  45. Cheb_vector & Cheb_vector::prime()
  46. {
  47.     Cheb_vector & deriv = *new Cheb_vector(size(),-1);
  48.     for (int n = 0; n < size(); n++) {
  49.     double summation = 0;
  50.     for (int p = n+1; p < size(); p++)
  51.         if ((p+n)%2)    /* if p+n is odd, we get a remainder */
  52.         summation += p * (*this)[p];
  53.     deriv[n] = 2.0*summation/C(n);
  54.     summation = 0;
  55.     }
  56.     return deriv;
  57. }
  58.  
  59.  
  60. /* 
  61. -*++ phys_vector::Chebyshev(): Transform to Chebyshev space
  62. ** 
  63. ** (*++ history: 
  64. **     17 Dec 87    Bruce Eckel    Creation date
  65. ** ++*)
  66. ** 
  67. ** (*++ detailed: 
  68. ** ++*)
  69. */
  70.  
  71. Cheb_vector & phys_vector::Chebyshev()
  72. {
  73.     Cheb_vector & cheb = *new Cheb_vector(size(),-1);
  74.     for (int mode = 0; mode < size(); mode++) {
  75.     double fact = M_PI * mode / size();
  76.     double sum = 0;
  77.     for ( int i = 0; i < size(); i++)
  78.         sum += (*this)[i] * cos( ((double)i + 0.5) * fact);
  79.     cheb[mode] = 2 * sum/size();
  80.     }
  81.     return cheb;
  82. }
  83.  
  84. /* 
  85. -*++ Cheb_vector::operator+(): function for Chebyshev vector addition
  86. ** 
  87. ** (*++ history: 
  88. **     15 Jan 88    Bruce Eckel    Creation date
  89. ** ++*)
  90. ** 
  91. ** (*++ detailed: 
  92. ** ++*)
  93. */
  94.  
  95. Cheb_vector & Cheb_vector::operator+(Cheb_vector & arg)
  96. {
  97.     if(size() != arg.size()) 
  98.     error("Cheb_vectors must be same size to add!");
  99.     Cheb_vector & sum = *new Cheb_vector(size(),0);
  100.     for(int i = 0; i < size(); i++)
  101.     sum[i] = (*this)[i] + arg[i];
  102.     return sum;
  103. }
  104.  
  105.  
  106. /* 
  107. -*++ Cheb_vector::operator-(): vector subtraction
  108. ** 
  109. ** (*++ history: 
  110. **     16 Jan 88    Bruce &    Creation date
  111. ** ++*)
  112. ** 
  113. ** (*++ detailed: 
  114. ** ++*)
  115. */
  116.  
  117. Cheb_vector & Cheb_vector::operator-(Cheb_vector & arg)
  118. {
  119.     if(size() != arg.size()) 
  120.     error("Cheb_vectors must be same size to subtract!");
  121.     Cheb_vector & sum = *new Cheb_vector(size(),0);
  122.     for(int i = 0; i < size(); i++)
  123.     sum[i] = (*this)[i] - arg[i];
  124.     return sum;
  125. }
  126.  
  127. /* 
  128. -*++ Cheb_vector::operator*(): multiply a Cheb_vector by a double
  129. ** 
  130. ** (*++ history: 
  131. **     16 Jan 88    Bruce &    Creation date
  132. ** ++*)
  133. ** 
  134. ** (*++ detailed: 
  135. ** ++*)
  136. */
  137.  
  138. Cheb_vector & Cheb_vector::operator*(double & arg)
  139. {
  140.     Cheb_vector & result = *new Cheb_vector(size(),0);
  141.     for(int i = 0; i < size(); i++)
  142.     result[i] = (*this)[i] * arg;
  143.     return result;
  144. }
  145.  
  146.