%%
% We create two vectors for data about the population of China. The first
% has the years of census data, the other has the numbers of millions of
% people.
year = (1980:10:2010)';  
population = [984.736; 1148.364; 1263.638; 1330.141];

%%
% It's convenient to measure time in years since 1980.
t = year - 1980;
y = population;

%%
% Now we have four data points $(t_1,y_1),\ldots,(t_4,y_4)$, so $n=4$ and
% we seek an interpolating cubic polynomial. We construct the associated
% Vandermonde matrix: 


%%
% To solve for the vector of polynomial coefficients, we use a backslash:


%%
% The algorithms used by the backslash operator are the main topic of next
% lectures. For now, observe that the coefficients of the cubic polynomial
% vary over several orders of magnitude, which is typical in this
% context. 

%% 
% Interpolating polynomial
f = @(c, t) c(1) + c(2)*t + c(3)*t.^2 + c(4)*t.^3

%%
% We can use the resulting polynomial to estimate the population of China
% in 2005:
estimate = f(c, 2005-1980);   % apply the 1980 time shift
disp(estimate)
%%
% The official figure is 1297.8, so our result is not bad. 

%%
% We can visualize the interpolation process. First, we 
% plot the data as points. We'll shift the t variable back to actual years. 
plot(1980+t, y, 'r*')

%%
% We want to superimpose a plot of the polynomial. In order to add to a
% plot, we must use the hold command:
hold on

%%
% To plot the interpolating polynomial, we create a vector with many points
% in the time interval using linspace.
tt = linspace(0, 30, 300)';   % 300 times from 1980 to 2010
yy = f(c, tt);          % evaluate the cubic
plot(1980+tt, yy)

%%
% Let's clear the figure (clf) and redo it, this time 
% continuing the curve outside of the original date range. We'll also 
% annotate the graph (using title, xlabel, ylabel, grid, and legend@ 
% to make its purpose clear.
clf   % clear figure
plot(1980+t, y, '.')
hold on
tt = linspace(-10, 50, 300)';   
plot(1980+tt, f(c,tt))
title('Population of China')
xlabel('year'); 
ylabel('population (millions)');
legend('data','interpolant','location','northwest');
grid on