% Script 3 for Week 11

% Comments on Chapra Problem 23.13
% The specified method ode45 is problematic here, and something is obviously 
% wrong with the solution, at least for the third function c3, of this 
% very stiff ode using ode45.
% We can fiddle with the ode45 options, seen below, to get different solution
% estimates. 

g=@(t,c) [-.013*c(1)-1000*c(1)*c(3); -2500*c(2)*c(3); -.013*c(1)-1000*c(1)*c(3)-2500*c(2)*c(3)]

% Notice that with initial conditions [c1(0),c2(0),c3(0)]=[1,1,0] the numerical
% estimate at time h (h is the step size used by the method for the first step) 
% for c3(h) will be negative and very close to -.013*h no matter how tiny h is. 
% After that the other two terms in the derivative of c3 (the third component 
% of g) will be positive so at the next step c3 will be less negative.
% The meaning of a negative c3 in this physical context is unclear.

% For some fairly substantial initial interval the solution to the third 
% component function c3 must behave like the solution to y' = -0.013 - 3500*y 
% which is a linear DE, which has exact solution

yexact=@(t) -0.013/3500 + 0.013/3500*exp(-3500*t);

% We should keep that in mind when we consider the validity of solutions. 
% c3 should look like the following, at least for a small initial interval

figure(26)
t=linspace(0,0.01);
y=yexact(t);
plot(t,y);

[u45,v45]=ode45(g,[0 50],[1 1 0]);  
figure(27)
plot(u45,v45(:,3))	
title('ode45 c3 coordinate')

[rk4,vrk4]=rk4sys(g,[0 50],[1 1 0], 0.001);  
figure(28)
plot(rk4,vrk4(:,3))	
title('rk4 c3 coordinate')

% rk4 does not seem to behave correctly. Initial c3 values seem positive. 
% ode45 just looks crazy.

[u45,v45]=ode45(g,[0 0.1],[1 1 0]);  
figure(29)
plot(u45,v45(:,3))	
title('ode45 c3 coordinate')

[rk4,vrk4]=rk4sys(g,[0 0.1],[1 1 0], 0.000001);  
figure(30)
plot(rk4,vrk4(:,3))	
title('rk4 c3 coordinate')

% The small t values are where the problem seems to lie. 

 OPTIONS = odeset('MaxStep',0.000001,'Stats', 'on')

[u45,v45]=ode45(g,[0 0.01],[1 1 0],OPTIONS);  
figure(31)
plot(u45,v45(:,3))	
title('ode45 with options c3 coordinate')


[rk4,vrk4]=rk4sys(g,[0 0.01],[1 1 0], 0.0000001);  
figure(32)
plot(rk4,vrk4(:,3))	
title('rk4 c3 coordinate')



%***********************************************************
% Extension of Chapra Problem 23.15

% The DE in problem 15 is 

% theta''=g/L*theta 

% and you should solve the problem with this DE.

% However ... the ACTUAL DE for this physical system on [0,pi] should be 

% theta''=g/L*sin(theta). (See the posted pdf for the derivation.)

% The given DE is the small-angle approximation, and we are not interested in 
% small angles only here.

% Following through .... rewrite the single 2nd order DE as the system

% V'=g/L*sin(T) 
% and T'=V

% with initial conditions 
% V(time 0)=0.25 and T(time 0)=0

%*****************************************
g=9.81
L=0.5
thetaAt0=0
thetaprimeAt0=0.25

f=@(t,c) [g/L*sin(c(2)); c(1)]

OPTIONS = odeset('MaxStep',0.001)
sol45=@(T) ode45(f,[0 T],[thetaprimeAt0 thetaAt0],OPTIONS);  
[u45,v45]=sol45(6);


figure(19)
plot(u45,v45(:,1))
title('time versus angular velocity (direct)')

xlabel('time, (s)')
ylabel('angular velocity (rad/sec)')

figure(20)
plot(u45,v45(:,2))
title('time versus angle (direct)')

xlabel('time, (s)')
ylabel('angle (rad)')

% Let's find the half-period of the motion, the time for angle pi.

label45=min(find(v45(:,2)>pi))
u45(label45(1))
v45(label45(1),:)


% produces output

label45 =        4475
ans =     1.1183
ans =    8.8624    3.1421

% So the model should switch from going faster to going slower a little 
% before time 1.1183 seconds.

% Running the code again twice more with a better interval and then finding a
% weighted average gives a very good estimate for the time when theta is pi. 
% I then plug the solution values on this interval into the function for 
% the solution to text problem 23.15.


[u45,v45]=sol45(1.12);
label45=min(find(v45(:,2)>pi));
TL=u45(label45(1))
TE=u45(label45(1)-1)
VL=v45(label45(1),:)
VE=v45(label45(1)-1,:)
C=(pi-VL(2))/(VE(2)-VL(2))
format long
Tpi= C*TE+(1-C)*TL
[u45,v45]=sol45(Tpi);
v45(length(v45),:)
format short
*******************************
% TL =    1.1183
% TE =    1.1181
% VL =    8.8624    3.1421
% VE =    8.8624    3.1399
% C =     0.2200
% Tpi =   1.118248809408893
% ans =   8.862420662550370   3.141592653546276
% (In finding C and Tpi, note the similarity to the shooting method of Ch. 24.)
% We have found the half-period of the motion, that produces angle pi 
% to 10 places. 
*******************************

% The following code mimics what you would HAVE to do to extend Chapra's 
% small angle approximation theta''=g/L*theta beyond angle pi. It uses 
% the symmetry of the solution required by conservation of energy to 
% "mirror" the first segment onto subsequent pi intervals. The sine 
% function takes care of that for us, but it works for us too!
% It uses the default options on ode45, as you would likely do in 
% the text problem.

[u45,v45]=sol45(Tpi);

N=length(v45)

t=zeros(1,6*N);
vel=zeros(1,6*N);
theta=zeros(1,6*N);

vel(1:N)= v45(:,1);
vel(2*N+1:3*N)=v45(:,1);
vel(4*N+1:5*N)=v45(:,1);

vel(N+1:2*N)=flip(v45(:,1));
vel(3*N+1:4*N)=flip(v45(:,1));
vel(5*N+1:6*N)=flip(v45(:,1));

theta(1:N)= v45(:,2);
theta(2*N+1:3*N)=2*pi+v45(:,2);
theta(4*N+1:5*N)=4*pi+v45(:,2);

theta(N+1:2*N)=2*pi-flip(v45(:,2));
theta(3*N+1:4*N)=4*pi-flip(v45(:,2));
theta(5*N+1:6*N)=6*pi-flip(v45(:,2));

for k=1:6
t(((k-1)*N+1) : (k*N))=(k-1)*Tpi+u45;
end

figure(33)
plot(t,vel)
title('time versus angular velocity (mirrored)')

xlabel('time, (s)')
ylabel('angular velocity (rad/sec)')

figure(34)
plot(t,theta)
title('time versus angle (mirrored)')

xlabel('time, (s)')
ylabel('angle (rad)')
*******************************

% This looks very similar to the solution you should get for Chapra's 
% problem 23.15  but with significantly longer period. 

%***********************************************************
Getting ode45 to equal the performance of rk4.m, and generally configuring options in the ode solvers.


odeset
          AbsTol: [ positive scalar or vector {1e-6} ]
          RelTol: [ positive scalar {1e-3} ]
     NormControl: [ on | {off} ]
     NonNegative: [ vector of integers ]
       OutputFcn: [ function_handle ]
       OutputSel: [ vector of integers ]
          Refine: [ positive integer ]
           Stats: [ on | {off} ]
     InitialStep: [ positive scalar ]
         MaxStep: [ positive scalar ]
             BDF: [ on | {off} ]
        MaxOrder: [ 1 | 2 | 3 | 4 | {5} ]
        Jacobian: [ matrix | function_handle ]
        JPattern: [ sparse matrix ]
      Vectorized: [ on | {off} ]
            Mass: [ matrix | function_handle ]
MStateDependence: [ none | {weak} | strong ]
       MvPattern: [ sparse matrix ]
    MassSingular: [ yes | no | {maybe} ]
    InitialSlope: [ vector ]
          Events: [ function_handle ]



OPTIONS = odeset(OLDOPTS,'NAME1',VALUE1,...)

Stats - Display computational cost statistics  [ on | {off} ]

MaxStep - Upper bound on step size  [ positive scalar ]
MaxStep defaults to one-tenth of the tspan interval in all solvers.

OPTIONS = odeset('MaxStep',0.001)

[TOUT,YOUT] = ode45(ODEFUN,TSPAN,Y0,OPTIONS)

f=@(t,x) [-2*x(1)+5*exp(-t) ; -x(1).*x(2).^2/2]

[a1,b1]=ode45(f,[0 .4],[2 4]);  

OPTIONS = odeset('MaxStep',0.1)
[a2,b2]=ode45(f,[0 .4],[2 4],OPTIONS);  


OPTIONS = odeset('MaxStep',0.01,'Stats', 'on')
[a3,b3]=ode45(f,[0 .4],[2 4],OPTIONS);  

OPTIONS = odeset('MaxStep',0.001,'Stats', 'on')
[a4,b4]=ode45(f,[0 .4],[2 4],OPTIONS);  

plot(a1,b1,a2,b2,a3,b3,a4,b4)
legend('default ode45','MaxStep 0.1','MaxStep 0.01','MaxStep 0.001')

% MaxStep 0.1 is terrible, but you have to zoom in quite a way 
% to distinguish the other three.

% More on Chapra Problem 24.2 using the shooting method.

% The shooting method solution using ode45 with default settings 
% reveals a lot of error relative to the analytic solution.
% We change some options below.

f2=@(t,x) [(x(2)) ;0.15* x(1)]

OPTIONS = odeset('MaxStep',0.001,'Stats', 'on')

[ x2,sol2]=ode45(f2,[0 10],[240 -93], OPTIONS);  
sol2(length(x2),:)
First= sol2(length(x2),2)

[ x2,sol2]=ode45(f2,[0 10],[240 -92],OPTIONS);  
sol2(length(x2),:)
Second= sol2(length(x2),2)

C=-Second/(First-Second)
G=C*(-93)+(1-C)*(-92)

[ x2,sol2]=ode45(f2,[0 10],[240 G],OPTIONS); 

sol2(length(x2),:)
sol2(1,:)
 
 A=-240*sinh(sqrt(0.15)*10)/cosh(sqrt(0.15)*10);
 B=240;
 F2an=@(x) A*sinh(sqrt(0.15)*x)+ B*cosh(sqrt(0.15)*x)
x2analytic=linspace(0,10);
y2analytic=F2an(x2analytic);

 
plot(x2analytic,y2analytic,x2, sol2(:,1))
title('position versus temperature')

xlabel('position on rod, (m)')
ylabel('temperature (K)')
legend('analytic','shooting with ode45')
legend('location','northwest')