home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / amiga / misc / icalc.lha / ICalc / Scripts / cheby.ic < prev    next >
Encoding:
Text File  |  1992-08-03  |  1.2 KB  |  63 lines
  1. # Routines for finding and evaluating Chebychev polynomials
  2. # with example of usage
  3. #
  4. # nb: routines assume array base is 1
  5. # nb: chebyfit takes a while...
  6.  
  7. NMAX = 40
  8. array ch_fnvals[1]    # array of function values for chebyfit (will be resized)
  9.  
  10. func chebyfit(a,b,@c,n,~f) = {
  11.     local bma, bpa, y, j, k, fac, sum
  12.     
  13.     resize(ch_fnvals, n)
  14.     bma = 0.5*(b-a)
  15.     bpa = 0.5*(b+a)
  16.  
  17.     for (k = 1; k <= n; k += 1)
  18.     {
  19.         y = cos(PI*(k-0.5)/n)
  20.         x = y*bma+bpa        # set up x
  21.         ch_fnvals[k] = f    # for function evaluation
  22.     }
  23.  
  24.     fac = 2/n
  25.     for (j = 1; j <= n; j += 1)
  26.     {
  27.         sum = 0
  28.         for (k = 1; k <= n; k += 1)
  29.             sum += ch_fnvals[k]*cos((PI*(j-1))*((k-0.5)/n))
  30.         c[j] = fac*sum
  31.     }
  32. }
  33.  
  34.  
  35. func chebyeval(a,b,@c,m,x) = {        # m is how many coeffs to use
  36.     local j, d, dd, y, y2, sv    # ...can (and should) be < NMAX
  37.  
  38.     if ((x-a)*(x-b) > 0) error(1)
  39.     d = 0
  40.     dd = 0
  41.     y = (2*x-a-b)/(b-a)
  42.     y2 = 2*y
  43.  
  44.     for (j = m; j >= 2; j -= 1)
  45.     {
  46.         sv = d
  47.         d = y2*d-dd+c[j]
  48.         dd = sv
  49.     }
  50.  
  51.     return(y*d-dd+0.5*c[1])
  52. }
  53.  
  54. # test functions out - find Chebychev poly for degree-arcsine
  55. array asc[NMAX]
  56. echo "Calculating Chebychev polynomial - be patient"
  57. chebyfit(-1,1,asc,NMAX,DEG*asin(x))
  58. echo "Evaluating some test values"
  59. chebyeval(-1,1,asc,20,-1)
  60. chebyeval(-1,1,asc,20,1)
  61. chebyeval(-1,1,asc,20,0.5)
  62. chebyeval(-1,1,asc,20,1/sqrt(2))
  63.