% This is the MATLAB diary from the lecture on the rail puzzle (NA review)
%
%
% the simple problem, if the rail bows up to make a tent with straight sides:
% (1+epsilon)^2 = 1^2 + d^2   <== solve for d
epsilon = 2/ 5280 % half mile units

epsilon =
  3.7879e-004

d = sqrt( 2*epsilon + epsilon^2)

d =
    0.0275

d * 5280/2

ans =
   72.6705

% the answer for the real rail problem should be between 0 & 72 feet.

% notation d = height
%          r = radius of circular arc
%        theta = angle of the arc
%          h = r - d = other leg of triangle with hypotenuse r  and angle theta

% We know the following equations
 % h = r - d
 % r theta = 5281/2 <-- let's make this our unit of measurement so that r theta = 1.
 % r sin theta = 5280/2 = 1 - epsilon
 % sin(theta) / theta = (1 - epsilon) / 1 % by eliminating r from previous two...
% Solve for theta:  sin(theta) = (1 - epsilon) * theta....

% fixed point iteration:
%   g(x) = sin(x) / (1 - epsilon) 
% theta is the solution to g(x) = x.

% Properties for fixed point iteration to converge:
%  g([a,b]) is contained in (a, b).
%  g(x) should be continuous on (a, b).
%      ===> there is fixed point in [a, b].

% |g'(x)| < 1 for all (a,b) then the FP is unique, and 
%             we can find the fixed point by iterating x = g(x).

x= 0.03

x =
    0.0300

format long

epsilon = 1/5281

epsilon =
    1.893580761219466e-004

x = sin(x) / (1 - epsilon)

x =
   0.03001372808339

x = sin(x) / (1 - epsilon)

x =
   0.03002059138874

x = sin(x) / (1 - epsilon)

x =
   0.03002745420136

x = sin(x) / (1 - epsilon)

x =
   0.03003431651987

x = sin(x) / (1 - epsilon)

x =
   0.03004117834289


% pretty slow; let's look at a spreadsheet for this iteration...
% Aitken's delta squared method helps some, and suggests that we try:

x = 0.0337

x = 
   0.033700000000

% The reason that we are converging so slowly is that the derivative is very close to one.

cos(x)/(1-epsilon)

ans =
   0.99980999208370

% Plotting this (which should usually be the first step) shows the same thing.

xx = 0:0.01:1;
plot(xx, sin(xx)/(1-epsilon))
hold on
plot(xx,xx,'r')


% remember: we want to find sin(theta)/theta = 1-epsilon....
%

% sin(x) = x - x^3/3!  + x^5/5! - x^7/7! + .....
% Taylor series are really useful.
% - they give polynomial approximations to hard functions
% - they give error bounds on how far the approximation is off...

%
% Taylor series f(x) at 0:
%  = f(0)/0! + f'(0)x/1! +f'2 (0) x^2 / 2! + ........

% sin(theta)/theta = 1 - theta^2/3! + theta^4/5! - theta^6/7! + ........
%
% sin(theta)/theta = 1-epsilon
%
% epsilon = theta^2/3! - theta^4 / 5! + theta^6 / 7! - ....
%
% lets look for theta^2 with the iteration
% theta^2 = 6 epsilon + (theta^2)^2/20 - (theta^2)^3/(7*6*5*4) 

theta2 = 0

theta2 =
     0

theta2 = 6* epsilon + (theta2)^2/20 - theta2^3/840

theta2 =
   0.00113640666810

theta2 = 6* epsilon + (theta2)^2/20 - theta2^3/840

theta2 =
   0.00113621302599

theta2 = 6* epsilon + (theta2)^2/20 - theta2^3/840

theta2 =
   0.00113621300399

theta2 = 6* epsilon + (theta2)^2/20 - theta2^3/840

theta2 =
   0.00113621300398

theta2 = 6* epsilon + (theta2)^2/20 - theta2^3/840

theta2 =
   0.00113621300398

% Much faster convergence; if we take the derivative we can find out why: it is about 1/10.

theta = sqrt(theta2)

theta =
   0.03370775880988

r = 1/theta

r =
  29.66676027440072

%%  Next :  solve d + h = r;  h^2 + (1-epsilon)^2 = r^2
%%
%%

% d = (r - h) (r+h)/(r+h) 
%   = ( r^2 - h^2 )/(r+h)  
%   = (1-epsilon)^2 / (r+h)

h = sqrt( r^2 - (1-epsilon)^2) 

h =
  29.64990799073509

d = (1-epsilon)^2/(r+h)

d =
   0.01685228366563

d * 5281/2 % back to units of feet

ans =
  44.49845501910026

% On thing I should have done differently is the evaluation of polynomials.
% 6*epsilon + x^2/20 - x^3/ 840 + x^4 / (9*8*840) <--- to evaluate a polynomial
% Horner's rule to evaluate polynomials
% (Synthetic division for dividing one polynomial by another)

(((x/9*8*840 +  -1/840)*x  + 1/20) *x + 0) *x + 6* epsilon   

ans =
   0.00215894062151


% You can also use Newton's method to determine theta as a zero of
%   f(x) = sin(x)/x - (1-epsilon) 
% 
% Newton's method is a fixed-point iteration with 
% newx = x - f(x)/f'(x)

% f'(x) = (x cos(x) - sin(x)) / x^2, so
% f(x)/f'(x) = (x sin(x) - x^2(1-epsilon)) / (x cos(x) - sin(x))
%
% I didn't use this in class because these subtractions may cause loss of precision, 
% since sin(x) ~=~ x and cos(x) ~=~ 1


x = 0.1

x =
   0.10000000000000

x = x - (x* sin(x) - x^2*(1-epsilon)) / (x *cos(x) - sin(x))

x =
   0.05566141357796

x = x - (x* sin(x) - x^2*(1-epsilon)) / (x *cos(x) - sin(x))

x =
   0.03803544681138

x = x - (x* sin(x) - x^2*(1-epsilon)) / (x *cos(x) - sin(x))

x =
   0.03395389846134

x = x - (x* sin(x) - x^2*(1-epsilon)) / (x *cos(x) - sin(x))

x =
   0.03370865076724

x = x - (x* sin(x) - x^2*(1-epsilon)) / (x *cos(x) - sin(x))

x =
   0.03370775882166

x = x - (x* sin(x) - x^2*(1-epsilon)) / (x *cos(x) - sin(x))

x =
   0.03370775880986

x = x - (x* sin(x) - x^2*(1-epsilon)) / (x *cos(x) - sin(x))

x =
   0.03370775880989

% It does, in fact, converge nicely.

diary off