% BRIEF SOLUTIONS TO NEEP271 PS4
%
% Fabian Waleffe 1998/11/04

%=========================================================
% PROBLEM 1
%
% Mass-Spring system was covered in class for general case
% Here all springs are identical so if x(i) is position of
% Mass number i and F(i) is the external force acting on it
% then balance of forces for that mass reads:
%
%  F(i) +K [x(i+1) - x(i) -L] - K [x(i) - x(i-1) - L] = 0
% or
%  K [-x(i-1) + 2 x(i) -x(i-1)] = F(i) = i (i-24) F
% with x(0)=0 and x(24)=24 L
%
% or  [-d(i-1) + 2 d(i) -d(i-1)] = i (i-24) F  with d(0)=d(24)=0
%
% and x(i) = i*L + d(i)
%
% L has dropped out because all the springs are identical!
% The d's are the DISPLACEMENTS from initial positions.
% The only relevant length scale in the problem is F/K
% so measure d in units of F/K : d -> F/K d
% 
% ==>   -d(i-1) +2 d(i) - d(i+1) = i(i-24)
% 

N=23; A=zeros(N); f=zeros(N,1);
for i=1:23,  A(i,i)=2;  f(i)=i*(i-24); end
for i=1:22,  A(i,i+1)=-1; A(i+1,i)=-1; end

d=A\f;

plot(d,'-o'), ylabel('K d/F'), xlabel('i')

%=========================================================
%  PROBLEM 2: see handout on pdf and random variables
%
% (1) prob to be in [a b]:
a=2; b=4; x0=1; s=3; 
prob=0.5*(erf((b-x0)/(sqrt(2)*s))-erf((a-x0)/(sqrt(2)*s)))

% (2) increasing samples

data=[];
for j=1:10,
   N=50*2^j; x=1+3*randn(N,1);
   data=[data;N, mean(x), std(x)];
end

Ns=data(:,1); means=data(:,2); stds=data(:,3);

loglog(Ns,abs(means-1),'-o',Ns,abs(stds-3)/3,'-o')

% The larger N the closer we are to mean=1 and std=3
% however everything is random, so we are only more likely
% to be close to those values when N is large.
% Statistical convergence to 1 and 3 is like sqrt(N),
% so it is slow on top of being only statistical

%=========================================================
% PROBLEM 3 was cancelled. See solutions on Web page.

%=========================================================
% PROBLEM 4: use Monte-Carlo integration
%
% Let's do 10 differents throwing-of-darts competitions:

vol=[]; nc=10;
for i=1:nc,
   N=1000; x=rand(N,1); y=-0.5+rand(N,1); z=rand(N,1);
   hem=x.^2+y.^2+z.^2; cyl=(x-0.5).^2+y.^2;
   vol=[vol; length(find(hem<1 & cyl<0.25))/N ];
end

[mean(vol), std(vol)]

exact=pi/3-4/9