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

clear 'all';

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

syms x y z
eq1 = exp(x^2+y^2-z^2)  == 9;
eq2 = x + y*x - z == 1;
eq3 = x+y+z  == 1
[x,y,z] = solve([eq1,eq2,eq3],[x,y,z])

% [2 POINTS]
%
%
% Problem 2: Plot the Pareto distribution
% P\left(x,\sigma ,\alpha \right)={\left(1+\frac{x}{\sigma }\right)}^{-\alpha }
%

clear 'all'
syms x sigma alpha real 
syms P(x,sigma, alpha)
P(x,sigma,alpha) = (1+x/sigma)^(-alpha);
figure
fplot(P(x,1,1.1),[0.0 3],'DisplayName','P(x,1,1.1)','Color',[0 0.447 0.741],...
     'LineWidth',2);
hold on
fplot( P(x,1,3),[0.0 3],'DisplayName','P(x,1,2)','Color',[0.85 0.325 0.098],...
     'LineWidth',2);
fplot(P(x,1,5),[0.0 3],'DisplayName','P(x,1,5)','Color',[0.929 0.694 0.125],...
     'LineWidth',2);
axes1 = gca;
xlabel({'x'},'FontAngle','italic','FontName','Times New Roman');
title('The Pareto distribution');
ylabel({'P(x,\sigma,\alpha)'},'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);

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

syms xi intvalue real
assume(sigma>0)
assume(alpha>1)
intvalue = simplify(int(P(xi,sigma,alpha),xi,[0,x]))

% [2 POINTS]
%
% Problem 4: Find the values of x that maximize Q\left(x\right)=\frac{\left(x^2 -1\right)}{x^4 +7}

syms Q(x) x xatmax
Q(x) = (x^2-1)/(x^4+7)
xatstationarypts = simplify(solve(diff(Q,x)==0,x))

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

Qatmax = Q(xatstationarypts)
% Depending on the version of matlab you are using the
% x values and Q values get stored in different places. Also some versions may
% return more than 3 roots for x.   The code below
% finds the index of the maximum value
[maxQ,index] = max(Qatmax)
double(Qatmax(index)) % The second and third ones are maxima...

% [2 POINTS]
%
%
%
%
% Problem 6: Solve the differential equation
% 6.1
% \frac{d\;V}{\textrm{dt}}=P+\eta V-\mu \left(S-V\right)

clear all
syms eta P S A t mu V0 t0 real
syms V(t)
IC = V(0)==V0
diffeq = diff(V(t),t) == P + eta*V(t) - mu*(S-V(t)) 
Vsol = dsolve([diffeq,IC],symvar(V(t)))

% [2 POINTS]
%
% 6.2 Solve for V_0

V0val = solve(subs(Vsol,t,t0)==0,V0)
Vsol = simplify(subs(Vsol,V0,V0val))
Vplot(t,mu) = subs(Vsol,[t0,eta,P,S],[10,0.04,800,1000000]);
figure
fplot(Vplot(t,0.002),[0 10],'DisplayName','\mu=0.002','Color',[0 0.447 0.741],...
     'LineWidth',2);
hold on
fplot(Vplot(t,0.02),[0 10],'DisplayName','\mu=0.02','Color',[0.85 0.325 0.098],...
     'LineWidth',2);
fplot(Vplot(t,0.034),[0 10],'DisplayName','\mu=0.034','Color',[0.929 0.694 0.125],...
     'LineWidth',2);
axes1 = gca;
xlabel({'Time (sec)'},'FontAngle','italic','FontName','Times New Roman');
title('Reserve for 10 year old term policy');
ylabel({'Reserve ($)'},'FontAngle','italic','FontName','Times New Roman');
set(axes1,'FontName','Times New Roman','FontSize',16,'XGrid','on','YGrid',...
    'on');
legend1 = legend(axes1,'show');
set(legend1,...
    'Position',[0.665983606557377 0.508333329529014 0.229166666666668 0.187022895474015],...
    'FontAngle','italic',...
    'FontSize',14);

% [2 POINTS + 1 POINT FOR LABELING PLOT]