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

clear 'all';

%
% Problem 1: Solve the equations
% \frac{1}{x}+\frac{1}{\left(y+z\right)}=-2\;\;\;\;\;\frac{1}{y}+\frac{1}{\left(x+z\right)}=-4\;\;\;\;\frac{1}{z}+\frac{1}{x}=-5\;\;\;\;\;\;\;\;

syms x y z
eq1 = 1/x + 1/(y+z)  == -2;
eq2 = 1/y + 1/(x+z) == -4;
eq3 = 1/z + 1/x  == -5;
[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 Gumbel distribution
% P\left(x,\mu ,\beta \right)=\frac{1}{\beta }\exp \left(-\left\lbrack \frac{\left(x-\mu \right)}{\beta }+\exp \left\lbrace -\frac{\left(x-\mu \right)}{\beta }\right\rbrace \right\rbrack \right)
%

clear 'all'
syms x xx beta mu real
syms P(x,mu, beta)
assume(beta>0); assume(mu>0);

P(x,mu,beta) = (1/beta)*exp( -( (x-mu)/beta+exp(-(x-mu)/beta) ) )
figure
fplot(P(x,0.5,1),[-5 15],'DisplayName','P(x,0.5,1)','Color',[0 0.447 0.741],...
     'LineWidth',2);
ylim([0, 0.5])
hold on
fplot( P(x,2,1),[-5 15],'DisplayName','P(x,2,1)','Color',[0.85 0.325 0.098],...
     'LineWidth',2);
fplot(P(x,2,4),[-5 15],'DisplayName','P(x,2,4)','Color',[0.929 0.694 0.125],...
     'LineWidth',2);
axes1 = gca;
xlabel({'x '},'FontAngle','italic','FontName','Times New Roman');
title('Gumbel distribution');
ylabel({'P(x,\mu,\beta)'},'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.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 \right)=\int_{-\infty }^x P\left(\xi ,\mu ,\beta \right)d\xi

syms xi real
syms S(x,mu, beta)
S(x,mu, beta) = simplify(int(P(xi,mu,beta),[-inf,x]))

% [1 POINTS]

figure
fplot(S(x,0.5,1),[0.0 15],'DisplayName','S(x,0.5,1)','Color',[0 0.447 0.741],...
     'LineWidth',2);
ylim([0, 1.2])
hold on
fplot( S(x,2,1),[0.0 15],'DisplayName','S(x,2,1)','Color',[0.85 0.325 0.098],...
     'LineWidth',2);
fplot(S(x,2,4),[0.0 15],'DisplayName','S(x,2,4)','Color',[0.929 0.694 0.125],...
     'LineWidth',2);
axes1 = gca;
xlabel({'x'},'FontAngle','italic','FontName','Times New Roman');
title('Cumulative Gumbel Distribution');
ylabel({'S(x,\mu,\beta)'},'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.6 0.2 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);

% [1 POINT + 1 POINT for labeling the plot]
%
% Problem 4: Find the value of x that maximizes P\left(x,\mu ,\beta \right)

xatmax = simplify(solve(diff(P(x,mu,beta),x)==0,x))  

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

Patmax = simplify(P(xatmax,mu,beta))

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

clear all
syms t mu lambda k N0 real
syms N(t)
IC = N(0)==N0
diffeq = diff(N(t),t) == mu*exp(-k*t)*N(t) - lambda*N(t)
Nsol(t,mu,k,lambda,N0) = simplify(dsolve([diffeq,IC],symvar(N(t))),'IgnoreAnalyticConstraints',true)

% [2 POINTS]
% 6.2

figure
fplot(real(Nsol(t,0.33,0.4,0.04,13.2)),[0 30],'DisplayName','\mu=0.33, k=0.4,\lambda=0.04','Color',[1 0 0],...
     'LineWidth',2);
hold on
fplot(real(Nsol(t,0.15,0.4,0.04,13.2)),[0 30],'DisplayName','\mu=0.15, k=0.4,\lambda=0.04','Color',[0 1 0],...
     'LineWidth',2);
fplot(real(Nsol(t,0.33,0.2,0.04,13.2)),[0 30],'DisplayName','\mu=0.33, k=0.2,\lambda=0.04','Color',[0 0 1],...
     'LineWidth',2);
fplot(real(Nsol(t,0.33,0.4,0.02,13.2)),[0 30],'DisplayName','\mu=0.33, k=0.4,\lambda=0.02','Color',[1 1 0],...
     'LineWidth',2);
axes1 = gca;
xlabel({'Time (days)'},'FontAngle','italic','FontName','Times New Roman');
title('Mammary Tissue Growth Model');
ylabel({'Cell mass (mg)'},'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]
%
%
%