%% Practical IV - Solutions to some problems

%% P1
% To get a random number between 2 and 10 to see what problem we want to
% attempt first, we call the randi method method to get a randon number 
% between 1 and 9, then add 1 to it (as we've already solved the first one):
firstProblem = randi(9)+1
%%
% Alternatively, to get a vector determining the problem order, we can find a
% random permuation of the vector 2:10 as follows:
problemOrder = randperm(9)+1


%% P2
% We create an anonymous function representing f as follows:
f = @(x) sin(x) + sin(x.^2);
%%
% We can then compute its integral as follows
fInt = integral(f,0,10)
%%
% To ensure we have obtained the value to the desired precision, we pass an
% extra tolerance argument to integral:
fInt12 = integral(f,0,10,'AbsTol',1e-12)
%%
% We check that the difference of the answers obtained is small
fInt- fInt12
%%
% to find the maximiser, we look for the minimiser of -f
x=fminbnd(@(x) -f(x),0,10);
fmax=f(x)
%%
% plot to check
fplot(f,[0,10])
hold on
plot(x,fmax,'r.','MarkerSize',10)
% see that this is only a local maximiser
%% repeat with smaller range of x
x2=fminbnd(@(x) -f(x),7.8,8.2);
fmax2=f(x2)
plot(x2,fmax2,'g.','MarkerSize',10)
grid
%%
% From the plot we see that we have roots around 2, 3, and 4, obtain them as follows:
r2 = fzero(f,2)
r3 = fzero(f,3)
r4 = fzero(f,4)
%%
% We confirm that those are indeed root by evaluating f at the points found
% and plotting
[f(r2) f(r3) f(r4)]
plot(r2,f(r2),'x',r3,f(r3),'x',r4,f(r4),'x')
hold off
%% P3
% Here, we need to combine integral, fzero, and the use of a double anonymous
% function. We need to create a function, I, which takes an input argument the
% parameter a, and returns the value of the integral
%   _
%  /
%  |  sin(x) + sin(ax^2) dx, from x = -1 to 1
% _/
%
% We achieve this as follows
I = @(a) integral(@(x) sin(x)+sin(a*x.^2),-1,1);
%%
% This can be interpreted as follows: For every value of a we pass to I,
% construct the anonymous function @(x) sin(x) + sin(a*x^2), and integrate it
% from -1 to 1.
%
% For example, here are some values of I
I(1)
I(4)
%%
% To find the value of a which gives I(a) = 1, we actually want to state the
% problem as a root-finding problem, that is, we want to use fzero to find a
% value a for which I(a) - 1 = 0. For that, we need to slightly redefine I as
% follows, subtracting 1 from the value returned
I = @(a) integral(@(x) sin(x)+sin(a*x.^2),-1,1)-1;
%%
% To get an idea where the root lies, we plot I as follows:
fplot(I,[0,5])
%%
% The roots seems to be close to a = 2, 3.2 so call fzero with that initial guess
a1 = fzero(I,2)
a2 = fzero(I,3.25)
%%
% Confirm that the values obtained are correct
I(a1)
I(a2)

%% P4
% MATLAB ships with U. S. census data for the national population every ten 
% years since 1790. Here we fit this with an exponential curve of the form 
%   c1 + c2 e^(c3*t)

%% Load the data

clc
clear
load census
whos

%% plot it

plot(cdate,pop,'-o')

%% Rescale time

t = (cdate-1790)/10;
plot(t,pop,'-o')

%% Make the fit

residual = @(c) norm( c(1) + c(2)*exp(c(3)*t) - pop ,2);
c = fminsearch(residual,[1 1 1])

%% Plot the fit

pp = c(1)+c(2)*exp(c(3)*t);
hold on
plot(t,pp,'--r','linewidth',2);
hold off

%% Extrapolate?

cdate2 = linspace(1790,2050,1000);
tt = (cdate2-1790)/10;
pp = c(1)+c(2)*exp(c(3)*tt);
figure
plot(cdate,pop,'-o'),  hold on
plot(cdate2,pp,'--r'); hold off

%% 500M?

f = @(t) c(1) + c(2)*exp(c(3)*t);
t_crit = fzero(@(t) f(t) - 500, 20)
y_crit = ceil(10*t_crit + 1790)


%% P5
% We call the roots command with the coefficients of the polynomial to find the
% roots of P as follows
coeffs = [1030 225 -1625 -290 610 50];
r = roots(coeffs)
%%
% We observe that we have five roots -- three real and two complex. 
%
% Using fzero, we create an anonymous function, and then try calling the method
% with different initial guesses.  Here, we can only obtain the real roots:
P = @(x) 50+610*x-290*x.^2-1625*x.^3+225*x.^4+1030*x.^5;
r1 = fzero(P,0)
r2 = fzero(P,.5)
r3 = fzero(P,1.1)

%% P7
% Here, we can use the 2D version of integral to find the answer. See >> help
% integral << for information about the calling sequence.
%
% Start by creating an anonymous function
f = @(x,y) y.*sin(x) + x.*cos(y);
fInt = integral2(f,pi,2*pi,0,pi)
% check - exact value is -pi^2
fInt+pi^2
%% P8
f=@(x)besselj(0,x);
fplot(f,[0,40]),grid
% 10th root is about 30
r10=fzero(f,30)
hold on, plot(r10,f(r10),'r.','MarkerSize',10), hold off

%% P9
% 
% To solve the IVP, create an anonymous function representing the differential
% equation
f = @(t,x)  t.*x;
%%
% Now call ode23, passing in the initial value of x, and the time interval we
% seek solutions on
[t, x] = ode23(f, [0 2], 1)

%%
% To find the value at an arbitrary time point, we can use the interp1 method,
% which interpolates the data and returns the value of the interpolant at the
% given point
x1 = interp1(t, x , 1)

%% repeat with ode45
[t, x] = ode45(f, [0,2], 1);
x2 = interp1(t, x, 1)

%% P10
% problem is equivalent to
% y1'=y2
% y2'=-y2*(1-y1^2)
dydx = @(x,y) [y(2); -y(2)*(1-y(1)^2)];
%%
% boundary conditions and initial guess
bcs = @(ya,yb) [ya(1)+1; yb(1)-1];
guess=@(x) [2*x-1, 2];
solinit = bvpinit(linspace(0,1,5),guess);
%% solve with bvp4c and plot
sol = bvp4c(dydx, bcs, solinit)
xx = linspace(0,1,501);
yy = deval(sol, xx);
plot(xx,yy,'LineWidth',2), legend('y','y''')

%% find where y=0.5
yc = @(c) interp1(xx,yy(1,:),c)-0.5;
xc = fzero(yc,0.6)

%% repeat with bvp5c
sol = bvp5c(dydx, bcs, solinit);
yy = deval(sol, xx);
yc = @(c) interp1(xx,yy(1,:),c)-0.5;
xc = fzero(yc,0.6)