%
% BRIEF SOLUTIONS TO PS5 
%
% Fabian Waleffe 1998/11/17
%
%===========================================================
% PROBLEM 1: SEE HANDOUT on "That crazy integral"
%            (SOlution of PS4 P3)
%            THINK POLAR!! area = integral r^2/2 dtheta
% compute from 0 to pi/2, then mult by 4
% answer is FULLY converged with N=18. WOW!

N=10; format long g, format compact
for i=1:6,
N=N+2;
t=[0:1/N:1]*pi/2; r=(cos(t).^4+sin(t).^4).^(-1/4); 
[N, trapz(t,r.^2/2)]
end

Area=4*trapz(t,r.^2/2)

%===========================================================
% PROBLEM 2: SEE HANDOUT on "That crazy integral"
%
% There is a singularity at x=0 but it is integrable.
%
% singularity is x^(-2/3)
% but x^(-2/3) dx = 3 d(x^(1/3)) so let u=x^(1/3) or x=u^3
% then x^(-2/3) dx = 3 du  and cos(x)=cos(u^3)
%

 N0=20;   format compact, format long g
 for i=0:3,
   N=N0*2^i;
   u=[0:1/N:1]*(pi/2)^(1/3);
   f=3*cos(u.^3);
   [N,trapz(u,f)]
 end

format 
integral=trapz(u,f)   %% Can trust first 4 digits

%===========================================================
% PROBLEM 3: SEE HANDOUT on Numerical quadrature
%
% You must have downloaded the "velo.dat" file as well
% as the "splineq.m" function and have told Matlab where
% to find them. If it is in the directory where Matlab 
% is running then 

load velo.dat, t=velo(:,1); v=velo(:,2);

distance=splineq(t,v)    
notsogooddistance=trapz(t,v)

%===========================================================
% PROBLEM 4: 
%
% rewrite as 1st order system: let y(1)=A, y(2)=A' then
% A'' + sin(A)=0  ------> dy/dt = [y(2); -sin(y(1))]

% The following function should be saved as file "FWpend.m" in you directory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dydt=FWpend(t,y)

% pendulum system A''+sin(A)=0 with y=[A; A']
% FWpend(t,y) returns dy/dt

dydt = [y(2); -sin(y(1))];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Let's integrate it from t=0 to t=8 with initial conditions A=pi/2, A'=0

A0=pi/2;
[t,y]=ode23('FWpend',[0 8],[A0; 0]);

% harmonic oscillator solution (see handout) is A0 cos(t)

Aharm=A0*cos(t);

% crazy A(t): 

sa0=sin(A0/2); m=sa0^2; d=ellipke(m);  Acrazy= 2*asin(sa0*ellipj(t+d,m));

% Let's plot them all on the same plot to compare:

plot(t,y(:,1), t, Aharm, t, Acrazy,':o'), xlabel('t'), ylabel('A'), title('pendulum') 

% that crazy function is the exact solution!! ode23 calculates it well.
% The harmonic oscillator is a good approx only for small initial A0.

% To determine period: look where y crosses zero, interpolate between t's. 
% Find index of 1st crossing using min(find(y<0)) then linear interpolation:
% y(t) ~ ym + slope * (t-tm) where slope = (yp-ym)/(tp-tm).

ip=min(find(y<0)); yp=y(ip); ym=y(ip-1); tp=t(ip); tm=t(ip-1);
tzero=tm - ym/(yp-ym)*(tp-tm);
approxperiod=4*tzero;  exactperiod=4*d; [approxperiod,exactperiod]

% You could recompute to try to get a better period, but one run
% followed by a linear interpolation gives a pretty good answer
% good enough for full credit. A0=pi/4 is same logic of course.