1 %{
2 Chebyshev Sine approximation
3 Octave 5.1.0
4 John Bryan 2019
5 %}
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=double(cos(0));
28 else
29 maxfunction=double(sin(pi/2));
30 endif
31 factor1=maxfunction;
32 factor2=double(1/factorial(n+1));
33 factor3=double(1/2^n);
34 factor4=double(((b-a)/2)^(n+1));
35 w=double(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 = double(a(1));
72 for i = 2:nplus1
73 result = double(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=double(zeros(n+1,n+1));
91 to_poly(1,n+1)=double(1);
92 for i=2:(n+1)
93 to_poly(i,n+2-i:n+1)=double(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=double(zeros(1,n+1));
111 % equation 2 implementation
112 x=double((b+a)/2+((b-a)/2)*cos((2*k+1)*pi/(2*n+2)));
113 % equation 3 implementation
114 c(1)=double((1/(n+1))*sum(f(x)));
115 % equation 4 implementation
116 factor=double(2/(n+1));
117 for j=1:n
118 summation=double(0);
119 for k=0:n
120 term = double(f(x(k+1)))*double(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 = double(0.0);
138 for i=1:(n+1)
139 chebpoly = chebpoly + double(c(i) * to_poly(i,:));
140 endfor
141 return;
142 endfunction
143
144
145
146 % f is a functionhandle for function being approximated
147 f=@sin;
148 % narray are the n values
149 % n=7 for single precision
150 % n=13 for double precision
151 narray=double([7,11,13]);
152 % a is domain lower limit
153 a=double(0);
154 % b is domain upper limit
155 b=double(pi/2);
156 % numtests is number of points in [a,b] domain tested.
157 numtests=double(129);
158 % x are the domain points tested.
159 x=double(pi*[0:(numtests-1)]/256);
160 for p=1:length(narray)
161 n=narray(p);
162 nplus1=n+1;
163 % to_poly are chebyshev polynomials
164 to_poly=double(calculatepolynomials(n));
165 % c are n+1 chebyshev coefficients
166 c=double(getcoefficients(f,n,a,b));
167 chebpoly=double(GetChebPoly(c,to_poly,n));
168 signeddifference=double(zeros(1,numtests));
169 t0=clock ();
170 for i=1:numtests
171 % y2 is built-in sine output
172 y2=double(sin(x(i)));
173 x2=x(i);
174 % parameter is change of variable due to using [a,b], not [-1,1]
175 parameter=double((2*(x2-a)/(b-a))-1);
176 % w is chebyshev approximation
177 w=NestedPolyEval(chebpoly, parameter, nplus1);
178 difference(i)=double(abs(double(w-y2)));
179 signeddifference(i)=double((w-y2));
180 endfor
181 elapsed_time=etime (clock (),t0)
182 maxdifference=double(max(difference))
183 TheoreticalAbsMaxError=double(TheoreticalMaxError(n,b,a))
184 PlotFunction(n,x,signeddifference,TheoreticalAbsMaxError,p)
185 endfor