1 %{
  2   Chebyshev Sine approximation.  
  3   Octave 5.1.0.
  4   John Bryan 2019.
  5   Retested in Octave 9.1.0 2024
  6  %}
  7 
  8 
  9  clear all;
 10  close all;
 11  warning('off','all');
 12  format long;
 13  pkg load symbolic;
 14  echo off;
 15 
 16 
 17  function w=TheoreticalMaxError(n,b,a)
 18     % Equation 6 implementation
 19     % Calculate theoretical max error bound for sine approx. on [0,pi/2].
 20     f = @(x) sin(x);
 21     % make anonymous function into a symbolic formula  
 22     syms x;
 23     ff = f(x);
 24     % calculate the (n+1) derivative of the function  
 25     ffdnp1 = diff(ff, x, n+1);
 26     if (isequal(abs(ffdnp1),abs(cos(x))))
 27        maxfunction=cos(0);
 28     else
 29        maxfunction=sin(pi/2);
 30     endif
 31     factor1=maxfunction;
 32     factor2=1/factorial(n+1);
 33     factor3=1/2^n;
 34     factor4=((b-a)/2)^(n+1);
 35     w=factor1*factor2*factor3*factor4;
 36     return;
 37  endfunction
 38 
 39 
 40  function PlotFunction(n,x,signeddifference,absmaxerror,p)
 41      % plots approximation error and error bound eq 6.
 42      % p is the figure number.
 43      % x is domain points tested.
 44      % signeddifference is approximation error
 45      % absmaxerror is eq 6 max error bound.
 46      figure (p);
 47      plot(x, signeddifference, 'b');
 48      hold on
 49      plot(x, absmaxerror,'r');
 50      hold on
 51      plot(x, -absmaxerror,'r');
 52      handle = legend (["n=" num2str(n) "approximation error"], ["n=" num2str(n) "theoretical max error"]);
 53      legend (handle, "location", "northeastoutside");
 54      legend boxoff;
 55      set(gcf,'InvertHardCopy','off');
 56      filename=["cheb_plot" int2str(p) ".png"];
 57      axis([0,pi/2]);
 58      print (gcf, "-S500,500", filename, '-dpng');
 59      return
 60  endfunction
 61 
 62 
 63  function result = NestedPolyEval(a, x, nplus1)
 64     % function performs nested polynomial evaluation;
 65     % input a=polynomial of order n;
 66     % a(1) is nth order term coefficient;
 67     % a(n+1) is zeroth order term;
 68     % input x is parameter to polynomial;
 69     % input nplus1 is n+1, the length of the polynomial;
 70     % output is result of nested polynomial evaluation;
 71     result = a(1);
 72     for i = 2:nplus1
 73       result = result*x + a(i);
 74     endfor
 75  endfunction
 76 
 77 
 78  function to_poly=calculatepolynomials(n)
 79     % equation 5 implementation.
 80     % calculate Chebyshev polynomials.
 81     % convert from symbols for speedup.
 82     syms x;
 83     syms to;
 84     to(1)=1;
 85     to(2)=x;
 86     for i=3:(n+1)
 87        to(i)=2*x*to(i-1)-to(i-2);
 88        to(i)=expand(to(i));
 89     endfor
 90     to_poly=zeros(n+1,n+1);
 91     to_poly(1,n+1)=1;
 92     for i=2:(n+1)
 93        to_poly(i,n+2-i:n+1)=sym2poly(to(i));
 94     endfor
 95     to_poly;
 96     return
 97  endfunction
 98 
 99 
100  function c=getcoefficients(f,n,a,b)
101    % implements equations 2, 3, 4.  
102    % returns Chebyshev coefficients c.
103    % f is functionhandle.
104    % n is order of Chebyshev polynomial.
105    % a is lower domain limit.
106    % b is upper domain limit.
107    % linspace(start,end,number of elements).
108    k=linspace(0,n,n+1);
109    % c is coefficients.
110    c=zeros(1,n+1);
111    % equation 2 implementation
112    x=(b+a)/2+((b-a)/2)*cos((2*k+1)*pi/(2*n+2));
113    % equation 3 implementation
114    c(1)=(1/(n+1))*sum(f(x));
115    % equation 4 implementation
116    factor=2/(n+1);
117    for j=1:n
118       summation=(0);
119       for k=0:n
120           term = f(x(k+1))*cos(j*(2*k+1)*pi/(2*n+2));
121           summation = summation + term;
122       endfor
123       summation=factor*summation;
124       c(j+1)=summation;
125    endfor
126    c;
127    return;
128  endfunction
129 
130 
131  function chebpoly = GetChebPoly(c,to_poly,n)
132  % equation 1, Chebyshev polynomial approximation equation.
133  % n=degree of polynomial.
134  % c=Chebyshev coefficients.
135  % to_poly are the n+1 Chebyshev polynomials.
136  % output chebpoly is the Chebyshev polynomial approximation.
137    chebpoly = 0.0;
138    for i=1:(n+1)
139       chebpoly = chebpoly + (c(i) * to_poly(i,:));
140    endfor
141    return;
142  endfunction
143 
144 
145 
146  function test()
147     % f is a functionhandle for function being approximated
148     f=@sin;
149     % narray are the n values
150     % n=7 for single precision
151     % n=13 for double precision
152     % n=12 results are a little noisy; a few results reach theoretical bound.
153     % n=13 results are a little more noisy; more results exceed theoretical bound.
154     % n=11 results are clean.
155     narray=[7,11];
156     % a is domain lower limit
157     a=(0);
158     % b is domain upper limit
159     b=(pi/2);
160     % numtests is number of points in [a,b] domain tested.
161     numtests=129;
162     % x are the domain points tested.
163     x=(pi*[0:(numtests-1)]/256);
164     for p=1:length(narray)
165       n=narray(p);
166       nplus1=n+1;
167       % to_poly are chebyshev polynomials
168       to_poly=calculatepolynomials(n);
169       % c are n+1 chebyshev coefficients
170       c=getcoefficients(f,n,a,b);
171       chebpoly=GetChebPoly(c,to_poly,n);
172       signeddifference=zeros(1,numtests);
173       t0=clock ();
174       for i=1:numtests
175           % y2 is built-in sine output
176           y2=sin(x(i));
177           x2=x(i);
178           % parameter is change of variable due to using [a,b], not [-1,1]
179           parameter=(2*(x2-a)/(b-a))-1;
180           % w is chebyshev approximation
181           w=NestedPolyEval(chebpoly, parameter, nplus1);
182           difference(i)=abs((w-y2));
183           signeddifference(i)=w-y2;
184       endfor
185       elapsed_time=etime(clock (),t0);
186       maxdifference=max(difference);
187       TheoreticalAbsMaxError=TheoreticalMaxError(n,b,a);
188       fprintf('\n n=%d Maximum difference: %1.15e \n',n,maxdifference);
189       fprintf('\n n=%d Theoretical absolute maximum error: %1.15e \n',n,TheoreticalAbsMaxError);
190       PlotFunction(n,x,signeddifference,TheoreticalAbsMaxError,p);
191     endfor
192   endfunction
193 
194 test()