%% ENGN0040 Dynamics and Vibrations: Solutions to HW1, Spring 2018
% [13 POINTS TOTAL]
%

clear 'all';  % Good idea to do this if we re-use variables in different problems

%
% Problem 1: Solve the equations
% x^2 +y^3 =4\text{      }x+y-x\text{ }y=1\text{         }

syms x y
eq1 = x^2+y^3==4;
eq2 = x + y - x*y==1;
[x,y] = solve([eq1,eq2],[x,y])

% [1 POINT]
%
%
% Problem 2: Plot the Weibull distribution
% P\left(x,\gamma \right)={\gamma x}^{\gamma -1} \exp \left(-x^{\gamma } \right)
%

clear 'x'
syms x gamma P
P = @(x,gamma) gamma*x^(gamma-1)*exp(-x^gamma);
fplot(P(x,0.5),[0.05 3],'DisplayName','P(x,0.5)','Color',[0 0.447 0.741],...
     'LineWidth',2);
hold on
fplot( P(x,2),[0.05 3],'DisplayName','P(x,2)','Color',[0.85 0.325 0.098],...
     'LineWidth',2);
fplot(P(x,5),[0.05 3],'DisplayName','P(x,5)','Color',[0.929 0.694 0.125],...
     'LineWidth',2);
axes1 = gca;
xlabel({'x'},'FontAngle','italic','FontName','Times New Roman');
title('The Weibull distribution function');
ylabel({'P(x,m)'},'FontAngle','italic','FontName','Times New Roman');
set(axes1,'FontName','Times New Roman','FontSize',20,'XGrid','on','YGrid',...
    'on');
legend1 = legend(axes1,'show');
set(legend1,...
    'Position',[0.665983606557377 0.708333329529014 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);

% [1 POINT + 1 POINT for labeling the plot]
%
% Problem 3:  Evaluate
% \int_0^x P\left(\xi ,\gamma \right)d\xi

syms xi intvalue
intvalue = int(P(xi,gamma),xi,[0,x])

% [1 POINT]
%
% Problem 4: Find the values of x that maximize P(x,5)

syms xatmax
xatmax = simplify(solve(diff(P(x,5),x)==0,x)) 
xatmax(2)  % We know this one is the max from the plot

% [1 POINT]
%
% Problem 5: Find the maximum value of P

Patmax = subs(P(x,5),x,xatmax(2))
double(Patmax)   % A numerical answer (not required)

% [1 POINT]
%
% Problem 6: Find the values of x that maximize P for arbitrary gamma.  (here we will look for local maxima not absolute maxima)
% There are many ways to solve this problem.   We can tell Matlab that we want solutions
% for real gamma, and only want real solutions, by using

assume(gamma,'real');     % We want real valued gamma
assume(xatmax,'real');        % We want real solutions
xatmax = solve(diff(P(x,gamma), x)==0, x) %

% This gives two solutions - Matlab tells us that the solutions are valid only for special values of \gamma.  The first solution is valid if \gamma >2 the second is valid for \gamma >1.   So we can clean things up a bit by doing the calculation again with these assumptions, and checking whether the solutions correspond to maxima, minima or inflection points by computing the second derivative.  Here is the first solution (zero)

clear 'gamma'
syms gamma
assume(gamma,'real')
assumeAlso(gamma>2)
xatmax = solve(diff(P(x,gamma), x)==0, x);
sol1 = xatmax(1)
secondderiv = simplify(subs(diff(P(x,gamma),x,2),x,sol1))

% This is positive (and infinite) or zero for \gamma >2 so it is a minimum or an inflection point and of no interest.
% Now the second

clear 'gamma'
syms gamma real
assume(gamma>1)
xatmax = solve(diff(P(x,gamma), x)==0, x);
sol2 = xatmax(2)
secondderiv = simplify(subs(diff(P(x,gamma),x,2),x,sol2))

% This gives a negative second derivative so its a maximum.   We conclude that {\left(1−\frac{1}{\gamma }\right)}^{\frac{1}{\gamma }} is the only solution, and it is valid for \gamma >1
%
% We can also solve the problem with the 'Returnconditions' variable in the solver

clear 'xatmax'
clear 'gamma'
syms gamma real
sols = solve(diff(P(x,gamma), x)==0, x,'ReturnConditions',true);
sols.x
sols.parameters
sols.conditions

% Matlab tells us what we just found out, but its a bit hard to follow the notation.  We see three solutions, the first of which is valid for \gamma >2, the second is valid for 1<\gamma <2, the third is valid for some weird combinations of conditions.  We need the integer parameter k=0 for the complex exponential in the last case to be real, so we can choose this one out of the family

solsnew = subs(sols.x,sols.parameters,0)
conditions = simplify(subs(sols.conditions,sols.parameters,0))

% This tells us that the first solution is valid for \gamma >2 - we could proceed as before to show that this root is a minimum or an inflection.   The second solution is valid for 1<\gamma <2, and the third (which is the same as the second) is valid for  \gamma >2  or  \gamma <0 .  We can check out the second derivative for \gamma >1 as before to conclude that {\left(1−\frac{1}{\gamma }\right)}^{\frac{1}{\gamma }} is indeed a maximum for \gamma >1 .  Solutions with  \gamma <0 are of no practical interest but we can check out the math and see whether the solution is min or max for this case

clear 'gamma'
syms gamma real
assume(gamma<0)
xatmax = solve(diff(P(x,gamma), x)==0, x)
secondderiv = simplify(subs(diff(P(x,gamma),x,2),x,xatmax))

% The second derivative is positive for  \gamma <0 and therefore a minimum and of no interest....
%
% Graders: there will be many solutions to this problem.   Please give 1.5 points for using the 'assume' approach and getting a correct answer, or using the 'returnconditions' and getting the right answer, and give 2 points for a more in-depth discussion along the lines presented here.  If solutions give absolute maxima that's OK too.
% [2 POINTS]
%
%
% Problem 7: Solve the differential equation
% \frac{\mathit{df}}{\mathit{dt}}=\mathit{rf}\left(t\right)-d\left(f\left(t\right)+\mathit{cf}{\left(t\right)}^2 \right)-H

syms diffeq f(t) t r d c H f0
f(t) = simplify(dsolve(diff(f(t),t) == r*f(t) - d*(f(t)+c*f(t)^2)-H,f(0)==f0,symvar(f(t))))

% [2 POINTS]
%
% 7.2 Plot f(t) This is the solution with r=2 (updated version of HW)

fnum(t,H) = real(subs(f(t),[r,d,c,f0],[2,5/100,1/100,1000]));
figure
fplot(fnum(t,0),[0 6],'DisplayName','H=0','Color',[0 0.447 0.741],...
     'LineWidth',2);
hold on
fplot(fnum(t,1200),[0 6],'DisplayName','H=1200','Color',[0.85 0.325 0.098],...
     'LineWidth',2);
fplot(fnum(t,1451),[0 6],'DisplayName','H=1451','Color',[0.929 0.694 0.125],...
     'LineWidth',2);
axes1 = gca;
xlabel({'Time (months)'},'FontAngle','italic','FontName','Times New Roman');
title('Gone Fishing....');
ylabel({'Fish Population'},'FontAngle','italic','FontName','Times New Roman');
set(axes1,'FontName','Times New Roman','FontSize',20,'XGrid','on','YGrid',...
    'on');
legend1 = legend(axes1,'show');
set(legend1,...
    'Position',[0.665983606557377 0.408333329529014 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);

%
% This is the solution with r=1 (original version of HW)

fnum(t,H) = real(subs(f(t),[r,d,c,f0],[1,5/100,1/100,1000]));
figure
fplot(fnum(t,0),[0 6],'DisplayName','H=0','Color',[0 0.447 0.741],...
     'LineWidth',2);
hold on
fplot(fnum(t,1200),[0 6],'DisplayName','H=1200','Color',[0.85 0.325 0.098],...
     'LineWidth',2);
fplot(fnum(t,1451),[0 6],'DisplayName','H=1451','Color',[0.929 0.694 0.125],...
     'LineWidth',2);
axes1 = gca;
ylim( [0 2000] )
xlabel({'Time (months)'},'FontAngle','italic','FontName','Times New Roman');
title('Gone Fishing....');
ylabel({'Fish Population'},'FontAngle','italic','FontName','Times New Roman');
set(axes1,'FontName','Times New Roman','FontSize',20,'XGrid','on','YGrid',...
    'on');
legend1 = legend(axes1,'show');
set(legend1,...
    'Position',[0.665983606557377 0.408333329529014 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);

%
% [2 POINTS + 1 POINT FOR LABELING PLOT]   (Solutions with r=1 or r=2 are both fine)