%% ENGN0040 Dynamics and Vibrations: Solutions to HW1, Spring 2021
% [16 POINTS TOTAL]
%

clear 'all';

%
% Problem 1: Solve the equations
% \exp \left({\left(x\right.}^2 +y^2 -z^2 \right)=9\;\;\;\;\;\;x+x\;y-z=1\;\;\;\;x+y+z=1\;\;\;\;\;\;\;\;

syms x y z
eq1 = exp((x/2)^2-(y/4)^2+(z/6)^2)  == 9;
eq2 = log(x + y - z*x) == 0;
eq3 = x+y+z  == 1;
[x,y,z] = solve([eq1,eq2,eq3],[x,y,z]);
x = simplify(x)
y = simplify(y)
z = simplify(z)

% [2 POINTS]
%
%
% Problem 2: Plot the Gompertz-Makeham distribution
% P\left(x,\alpha ,\beta ,\lambda \right)=\left(\alpha e^{\beta x} +\lambda \right)\exp \left(-\lambda x-\frac{\alpha }{\beta }\left\lbrace e^{\beta x} -1\right\rbrace \right)
%

clear 'all'
syms x xx alpha beta lambda real 
syms P(x,alpha, beta, lambda)

P(x,alpha,beta,lambda) = (alpha*exp(beta*x)+lambda)*exp( -lambda*x - (alpha/beta)*(exp(beta*x)-1));
figure
fplot(P(x,5.e-05,0.085,0.00),[0.0 100],'DisplayName','P(x,5.e-05,0.085,0)','Color',[0 0.447 0.741],...
     'LineWidth',2);
ylim([0, 0.045])
hold on
fplot( P(x,5.e-04,0.085,0.00),[0.0 100],'DisplayName','P(x,5.e-04,0.085,0)','Color',[0.85 0.325 0.098],...
     'LineWidth',2);
fplot(P(x,5.e-05,0.085,0.001),[0.0 100],'DisplayName','P(x,5.e-05,0.085,0.001)','Color',[0.929 0.694 0.125],...
     'LineWidth',2);
axes1 = gca;
xlabel({'x (age)'},'FontAngle','italic','FontName','Times New Roman');
title('Gompertz-Makeham distribution');
ylabel({'P(x,\alpha,\beta,\lambda)'},'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.2665983606557377 0.708333329529014 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);

% [3 POINTS + 1 POINT for labeling the plot]
%
% Problem 3:  Evaluate
% S\left(x,\alpha ,\beta ,\lambda \right)=1-\int_0^x P\left(x,\alpha ,\beta ,\lambda \right)d\xi

syms xi 
syms S(x,alpha, beta, lambda)
S(x,alpha, beta, lambda) = 1-simplify(int(P(xi,alpha,beta,lambda),xi,[0,x]))

% [2 POINTS]

figure
fplot(S(x,5.e-05,0.085,0.00),[0.0 100],'DisplayName','S(x,5.e-05,0.085,0)','Color',[0 0.447 0.741],...
     'LineWidth',2);
ylim([0, 1.2])
hold on
fplot( S(x,5.e-04,0.085,0.00),[0.0 100],'DisplayName','S(x,5.e-04,0.085,0)','Color',[0.85 0.325 0.098],...
     'LineWidth',2);
fplot(S(x,5.e-05,0.085,0.001),[0.0 100],'DisplayName','S(x,5.e-05,0.085,0.001)','Color',[0.929 0.694 0.125],...
     'LineWidth',2);
axes1 = gca;
xlabel({'x (age)'},'FontAngle','italic','FontName','Times New Roman');
title('Survival probability');
ylabel({'S(x,\alpha,\beta,\lambda)'},'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.2665983606557377 0.2 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);

% [2 POINTS + 1 POINT for labeling the plot]
%
% Problem 4: Find the values of x that maximize P\left(x,\alpha ,\beta ,\lambda \right)

xatstationarypts = simplify(solve(diff(P(x,alpha,beta,lambda),x)==0,x))

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

Patmax = simplify(P(xatstationarypts,alpha,beta,lambda))

% [1 POINT]
%
%
%
%
% Problem 6: Solve the differential equation
% 6.1
% \frac{d\;I}{\textrm{dt}}=\left(\lambda -\mu \right)I-\lambda I^2

clear all
syms t lambda mu I0
syms I(t)
IC = I(0)==I0
diffeq = diff(I(t),t) == (lambda-mu)*I(t) - lambda*I(t)^2 
Isol = simplify(dsolve([diffeq,IC],symvar(I(t))))

% [2 POINTS]
% 6.2

Iplot(t,mu) = subs(Isol,[I0,lambda],[0.05,1])
figure
fplot(Iplot(t,0.85),[0 30],'DisplayName','\mu=0.85','Color',[0 0.447 0.741],...
     'LineWidth',2);
 hold on
 fplot(Iplot(t,0.5),[0 30],'DisplayName','\mu=0.5','Color',[0.85 0.325 0.098],...
     'LineWidth',2);
 fplot(Iplot(t,0.1),[0 30],'DisplayName','\mu=0.1','Color',[0.929 0.694 0.125],...
     'LineWidth',2);
 axes1 = gca;
xlabel({'Time (days)'},'FontAngle','italic','FontName','Times New Roman');
title('SI Infection Model');
ylabel({'% Infected 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.6665983606557377 0.7 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);

% [2 POINTS + 1 POINT FOR LABELING PLOT]
%
%
%