1
2
3
4
5
6
7
8
9 clc;
10 echo off;
11 clear all;
12 close all;
13 warning off all;
14 format long;
15 pkg load symbolic;
16
17
18 function extremalist=Extremas(fun,a,b)
19 extremalist=[];
20 numpoints=200;
21 points=linspace(a,b,numpoints);
22 leftpoint=points(1);
23 rightpoint=points(2);
24 for (i=1:numpoints-1)
25 leftpoint=points(i);
26 rightpoint=points(i+1);
27 if ((fun(leftpoint)*fun(rightpoint))<0)
28 extremapoint=fzero(fun,[leftpoint,rightpoint]);
29 extremalist =[extremalist extremapoint];
30 endif
31 endfor
32 endfunction
33
34
35 function point=GetExtrema(remezmatrix,p3,a,b,f)
36 poly0=poly2sym(p3);
37 f=@(x) f(x)-poly0;
38 syms x;
39 ff=formula(f(x));
40 ffd=diff(ff);
41 df=function_handle(ffd);
42 extremalist=Extremas(df,a,b);
43 point=[a extremalist b];
44 endfunction
45
46
47 function p3=GetPolynomial(remezmatrix,point,n,f);
48 for i=1:(n+2)
49 for j=1:(n+1)
50 remezmatrix(i,j)=((point(i))^(j-1));
51 endfor
52 remezmatrix(i,n+2)=((-1)^(i-1));
53 endfor
54 fvector=(f(point))';
55 coeffperror=(remezmatrix\fvector);
56 p2=coeffperror(1:(n+1));
57 p3=fliplr(p2');
58 endfunction
59
60
61 function PlotFunction(xv,p3,point,count,f,g,count2)
62 xverror=f(xv)-polyval(p3,xv);
63 np2error=f(point)-polyval(p3,point);
64 figure;
65 plot(xv,xverror);
66 hold on;
67 line(xlim, [0,0], 'Color', 'k', 'LineWidth', 1);
68 hold on;
69 plot(point,np2error,"or");
70 handle = legend (["iteration " num2str(count) " error"]);
71 legend (handle, "location", "northeastoutside");
72 legend boxoff;
73 set(gcf,'InvertHardCopy','off');
74 title(g);
75 filename=["remez_plot" int2str(count2) int2str(count) ".png"];
76 print (gcf, "-S500,500", filename, '-dpng');
77 endfunction
78
79
80 function p3=RemezExchange(remezmatrix,p3,a,b,n,xv,unitythreshold,np2,f,g,count2)
81 count=1;
82 point=GetExtrema(remezmatrix,p3,a,b,f);
83 unityproximity=GetUnityProximity(point,np2,p3,f);
84 while (unityproximity>unitythreshold)
85 count=count+1;
86 p3=GetPolynomial(remezmatrix,point,n,f);
87 point=GetExtrema(remezmatrix,p3,a,b,f);
88 unityproximity=GetUnityProximity(point,np2,p3,f);
89 PlotFunction(xv,p3,point,count,f,g,count2);
90 endwhile
91 endfunction
92
93
94 function unityproximity=GetUnityProximity(point,np2,p3,f)
95 extremaerror=f(point)-polyval(p3,point);
96 smallest=abs(extremaerror(1));
97 largest=abs(extremaerror(1));
98 for i=2:np2
99 if (abs(extremaerror(i))<smallest)
100 smallest=abs(extremaerror(i));
101 elseif (abs(extremaerror(i))>largest)
102 largest=abs(extremaerror(i));
103 endif
104 endfor
105 unityproximity=largest/smallest;
106 endfunction
107
108
109 function test()
110 n=([4,8]);
111 np2=n+2;
112 a=([0,-2]);
113 b=([2,2]);
114 numtests=(400);
115 f={@(x)cos(exp(x)),@(x)sin(exp(x))};
116 fstring={"cos(exp(x)) N=4 approximation error";"sin(exp(x)) N=8 approximation error"};
117 unitythreshold=1.000005;
118 count=1;
119 count2=1;
120 for j=1:length(f)
121 remezmatrix=zeros(n(j)+2,n(j)+2);
122 iteration=1;
123 xv=linspace(a(j),b(j),numtests);
124 npoint=linspace(0,n(j)+1,n(j)+2);
125 remezmatrix=zeros(n(j)+2,n(j)+2);
126 point=0.5*(a(j)+b(j))+0.5*(b(j)-a(j))*cos(npoint*pi/(n(j)+1));
127 point=fliplr(point);
128 g=f{j};
129 gstring=fstring{j};
130 p3=GetPolynomial(remezmatrix,point,n(j),g);
131 PlotFunction(xv,p3,point,count,g,gstring,count2);
132 p4=RemezExchange(remezmatrix,p3,a(j),b(j),n(j),xv,unitythreshold,np2(j),g,gstring,count2);
133 count2=count2+1;
134 count=1;
135 endfor
136 endfunction
137
138
139 test()