%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solution to Problem 4:
% Sensitivity Analysis:
% When a Little Means a Lot
%
% Solve an ordinary differential equation  y' = a y, y(0) = 1
% over the interval 0 < t < 30 with  a =.006 (births and deaths)
% and a migration effect of 0 or bup=.003. 
%
% Then see how sensitive the equation is  by defining a year-by-year rate
% in the array myrate, where each entry is normally distributed with mean
% a and standard deviation tau.
%
% Repeat myrepetitions times and plot the results.
%
% exode.m     Dianne O'Leary 05/2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global myrate 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a = .006;
nrepetitions = 50;
length = 50;
tau = a/6;
bup = .003;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the problems with rate a and rate a+bup.
% We use ode45 with function values produced by myode.m,
% although it is easier to use the analytic solution.
% The global array myrate communicates the yearly growth rates
% to myode.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

myrate = a * ones(length+1,1);
[t0,y0]=ode45(@myode,[0:length],1);
myrate = (a+bup) * ones(length+1,1);
[t1,y1]=ode45(@myode,[0:length],1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Introduce yearly perturbations in the growth rate and
% solve the ode again (nrepetitions times). 
% Plot the results.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:nrepetitions,
   myrate = a + tau*randn(length+1,1) + bup*rand(length+1,1);
   [t,y]=ode45(@myode,[0:length],1);
   plot(t,y,'b')
   hold on
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Repeat, but order the high rates first.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:nrepetitions,
   myrate = a + tau*randn(length+1,1) + bup*rand(length+1,1);
   myrate = -sort(-myrate);
   [t,y]=ode45(@myode,[0:length],1);
   plot(t,y,'r')
   hold on
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Repeat, but order the low rates first.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:nrepetitions,
   myrate = a + tau*randn(length+1,1) + bup*rand(length+1,1);
   myrate = sort(myrate);
   [t,y]=ode45(@myode,[0:length],1);
   plot(t,y,'g')
   hold on
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the unperturbed results. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plot(t0,y0,'k')
plot(t1,y1,'k')
disp(sprintf('lower solution: %f %f',t0(end),y0(end)))
disp(sprintf('upper solution: %f %f',t1(end),y1(end)))

hold off
title('Solutions to the differential equation')
xlabel('t')
ylabel('y')
print -depsc ode_ex.eps