%% Lab 6 -- Leibniz's Wheel
%% Taylor polynomials
% In the Pre-Lab, we discuss using the Taylor Series for \ln(1 + x) given below
% \ln(1+x) = \sum_{n = 1}^\infty (-1)^{n+1} \frac{x^n}{n} = x - \frac{x^2}{2} + \frac{x^3}{3} - \frac{x^4}{4} + \ldots
% which has radius of convergence 1.
% To visualize approximating \ln(1 + x) with just a few terms of this series, the code below plots the natural log function (which remember in MATLAB is log and not ln) against its Maclaurin polynomial of degree N. We see that on the interval [-0.9,0.9] even the 5th or 6th degree Maclaurin polynomial gets you a pretty good approximation!

% Degree of the Maclaurin polynomial
N = 2

% The natural log function
f = @(x) log(1+x)

% The Maclaurin Polynomial of f
% First we define the array of terms that appear
g = @(x,n) (-1).^(n+1) .* x.^n ./ n
% Then we add it up to degree N
Tf = @(x) sum(g(x,transpose(1:N)))

fplot(f, [-0.9,0.9])
hold on;
fplot(Tf,[-0.9,0.9])
hold off;

%% Exercise 1
% We now finish the calculation from the Pre-Lab, where we found that the difference between the fifth degree Maclaurin polynomial for \ln(x + 1) and the full Maclaurin series above is at most
% \sum_{n = 6}^\infty |x|^n
% Now, this quantity is a geometric series, and so you can sum it. What is its sum? (Hint: \sum_{n = 6}^\infty |x |^n = \sum_{n = 0}^\infty |x|^6 \cdot |x|^n. This should be a power series you recognize.)
% What is the sum's maximum value on [-0.9,0.9]? What happens if you started with the 200th degree Maclaurin polynomial instead?
% (You may use a calculator or the space below to evaluate the numerical value.)



%
%% The natural log of 2
% Let's now try to use the Maclaurin polynomial to compute \ln(2).
% First, the value of \ln(2) is known to be

log(2)

% Let's see how well our Maclaurin series performs. First, \ln(2) = \ln(1+1) so we want to evaluate our Maclaurin polynomial at x = 1.
% Below we generate the sequence of numbers a_k, corresponding to the value of the kth degree Maclaurin polynomial of \ln(1+x) evaluated at x = 1. In other words, it is the sum of the first k terms of \sum (-1)^{n+1} x^n / n when x = 1.

% First we write down the sequence corresponding to the 
% Maclaurin series, with x = 1
fake_ln_term = @(n) (-1).^(n+1) ./ n
% The cumsum function takes cummulative sums, so the kth entry 
% in the array a corresponds to the value of the kth degree
% Maclaurin polynomnial evaluated at x = 1
a = cumsum(fake_ln_term(1:250))

% Let's plot the difference between the Maclaurin polynomial 
% approximation versus the real value, the x axis is the degree
% of the polynomial, the y axis is the difference

plot(a - log(2), 'ro')

% Let's zoom in near 200 to see what happens there
plot(175:225, a(175:225)- log(2), 'ro')

% Huh?! These differences seem much, much worse than what you saw in Exercise 1. (Here the differences is around 10^{-3}, compared to the much smaller value you found in Exercise 1 for the 200th degree polynomial.) (We will think about why later.)
%% Exercise 2
% By trial and error, estimate the smallest N such that the Nth degree Maclaurin polynomial Tf for \ln(1+x), when evaluated at x = 1, differs from \ln(2) by no more than 10^{-8}, the precision of most hand-held calculators. (You just need to find N within an order of magnitude.)

% The following computes the 153rd degree Maclaurin 
% polynomial for ln(1+x), evaluated at x = 1
value = sum(fake_ln_term(1:153))

% The difference between that and ln(2)

value - log(2)

%
%% The natural log of 2, take 2
% It turns out there's a different way of computing the natural log of 2, which requires fewer terms.
% Using the rule of logarithms, we remember that
% \ln(2) = - \ln(1/2)
% So we can also compute \ln(2) by computing \ln(0.5) and taking its negative.
% The following code implements this by applying the Maclaurin series formula for \ln(0.5) = \ln(1 + (-0.5)) = \sum (-1)^{n+1} (-0.5)^n / n and calculating up to degree N for various N.

% First we write down the sequence corresponding to the 
% Maclaurin series, with x = -0.5
fake_ln_term2 = @(n) (-1).^(n+1) .* (-0.5).^n ./ n
% The cumsum function takes cummulative sums, so the kth entry 
% in the array a corresponds to the value of the kth degree
% Maclaurin polynomnial evaluated at x = -0.5
b = cumsum(fake_ln_term2(1:250))

% Remember, b should be approximating the negative of ln(2)

% Let's plot the difference between the Maclaurin polynomial 
% approximation versus the real value

plot(-b - log(2), 'ro')

% Let's compare what happens between this method, against 
% the method used before exercise 2
% Here the plot has on the x axis the degree (up to 20) 
% The y axis shows the differences of the approximation to 
% ln(2), the method using ln(1+(-0.5)) in red, and the method 
% using ln(1+1) in blue

plot(-b(1:20)- log(2), 'ro')
hold on;
plot(a(1:20) - log(2), 'bo')
hold off;

% Looking at the plot, it looks like the approximation using \ln(1+(-0.5)) is approaching the real value much quicker (the approximation error being much smaller for the same degree of polynomial used).
%% Exercise 3
% Find the smallest degree N for which the Nth degree Maclaurin polynomial for f(x) = \ln(1+x) when evaluated at x = -0.5 differs from \ln(0.5) by less than 10^{-8. (For this one you should be able to find the exact value N with no problem.)

% The following computes the 53rd degree Maclaurin 
% polynomial for ln(1+x), evaluated at x = -0.5
value2 = sum(fake_ln_term2(1:53))

% The difference between that and ln(0.5)

value2 - log(0.5)

%
%% Comparison
% The N you found for Exercises 2 and 3 are roughly equal to how many basic computer operations that a calculator has to carry out in order to evaluate \ln(2) to 8 digit precision, using in Exercise 2 the approximation
% \ln(2) = \ln(1+1) \approx \sum_{n = 1}^N \frac{(-1)^{n+1}}{n}
% and in Exercise 3 the approximation
% \ln(2) = - \ln(0.5) = - \ln(1 + (-0.5)) \approx (-1)\cdot \sum_{n = 1}^N \frac{(-1)^{n+1} (-0.5)^n}{n}
% We see that the method in Exercise 3 works much better in practice, in the sense of needing much fewer operations before getting a good answer!
%% Exercise 4
% Discuss with your team why the method in Exercise 3 works much better than the method in Exercise 2.
% (Hint: it may help to look at the following plot of \ln(1+x) plotted against the 10th order Maclaurin polynomial. Why is the point x = 1 special?)

% Define f(x) = ln(1+x)
f = @(x) log(1+x)

% Define the 10th degree Maclaurin series of f
f_terms = @(x,n) (-1).^(n+1) .* x.^n ./ n
% Then we add it up to degree 10
Tf = @(x) sum(f_terms(x,transpose(1:10)))

fplot(f, [-0.6,1.2])
hold on;
fplot(Tf,[-0.6,1.2])
hold off;

%
%% Historical Remarks
% It may be surprising to know that computer science as a discipline traces back a history as long as calculus, considering that the modern electronic computers weren't invented until the 1930s, while calculus traces back to the 17th century.
% The title of this lab refers to an invention of Gottfried Wilhelm von Leibniz, who developed calculus contemporaneously and independently of Isaac Newton. Leibniz's wheel was an early design of a mechanical calculator.
% Leibniz, however, was way more ambitious than simply building calculators. He conceived of and designed a universal computing framework (the calculus ratiocinator) together with a formal language (the characteristica universalis), anticipating by three centuries modern computers and programming languages.
%% Modern Implementations
% In modern computer languages, the natural logarithm function is computed using a slightly modified version of the Taylor polynomial for \ln(x).
% In the C programming language as implemented in fdlibm library, for example, a 14th order polynomial is used: http://www.netlib.org/fdlibm/e_log.c
% In the Julia programming language, the implementation can be found here: https://github.com/JuliaLang/julia/blob/master/base/special/log.jl
% The modification made is for the purpose of "averaging out" the error. Recall from the course that Taylor's remainder theorem says:
%     If on the interval (-a,a) the k+1st derivative f^{(k+1)}(x) has absolute value no more than some number M, than the kth order Taylor polynomial approximates the function f by better than
% |T_kf(x) - f(x)| \leq \frac{M}{(k+1)!} x^{k+1}
% This remainder term is non-uniform over (-a,a): for x closer to zero the error is smaller. The modification made to the Taylor polynomial fine-tunes it by making the error near zero slightly bigger, while making the error slightly smaller for larger x.