%% Polynomial Data Fitting & Trigonometric Functions
%
% Given a set of points sampled form a trig function in R^2, find the polynomial that contains them.
%% Define the trig function you want to integrate
Given_f = @(x) sin(x);
%Given_f = @(x) exp(x.^2);
%Given_f = @(x) sin(pi*x.^2)/2;
%% We choose random points on the graph of $$\sin{\pi x^2 / 2}$.
k = 10; % number of points
x_given = -2+4*rand(1,k);
y_given = Given_f(x_given);
n = length(x_given);
%% The "Coefficient" Matrix
% The polynomial is a_0 + a_1*x + ... + a_m*x^m where m = n-1;
% This general form determines one equation for every given pair
M = ones(n,n+1);
for point = 1:n
for i = 1:n;
M(point,i) = x_given(point)^(i-1);
end
M(point,n+1) = y_given(point);
end
%% Solve the system
A = rref(M);
a = A(:,end);
%% The polynomial
% For example, y = a(1) + a(2)*x + a(3)*x.^2;
r = max(abs(x_given))+6; % How wide is the range of the plotting window
x = -r:.0001:r; % a sample of x values in the chosen range
s = [num2str(a(1))]; % make a string out of the number held in a(1)
% term by term, make a polynomial whose coeffs are taken from a
for i = 2:length(a);
% s is a character string that looks like the poly we want
s = [s,' + ',num2str(a(i)),'*x.^',num2str(i-1)];
end
% now we turn the character string into a set of y-values of the polynomial
disp('The approximating polynomial is: ')
disp(s)
PolyApprox = @(x) eval(s);
%% Plot the points and the curve in a psudo-animation
close all; % close all other figures
figure; % create a new figure
ry = max(abs(Given_f(x_given)))+20;
ylim([-ry ry]); % the max/min of the y-axis
xlim([-r r]); % the max/min of the x-axis
hold on % put all the plots onto the same figure
title('A trig function and its approximation'); % puts the string s into the title
% Plot the given points as large blue stars,
plot(x_given,y_given,'*b','MarkerSize',12);
disp('The given points are in blue')
legend('Given Points') % create a legend
pause
xx = [-r:.001:r]; % sample x-values, the domain of the function
yy = Given_f(xx);
plot(xx,yy,'sk','MarkerSize',1) % plot the exact function in black
disp('The given function is in black')
legend('Given Points', 'Given Function')
pause % pause the printing for greater appreciation
% Plot the graph of the polynomial along with the given points
plot(x,PolyApprox(x),'-r','MarkerSize',2);
plot(x_given,Given_f(x_given),'*b','MarkerSize',12);
disp('The polynomial approximation is in red')
legend('Given Points','Given Function', 'Polynomial Approx')
%% One could evaluate the exact and approximate integrals by calling the script CompareIntegrals afterward.