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

clear 'all';

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

syms x y z
eq1 = x^2+y^2 + z^2  == 4;
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 log-Cauchy distribution
% P\left(x,\sigma \right)=\left(\frac{1}{\pi x}\right)\left(\frac{\sigma }{{\mathrm{log}\left(x\right)}^2 +\sigma^2 }\right)
%

clear 'x'
syms x sigma real 
syms P(x,sigma)
P(x,sigma) = 1/(pi*x)* sigma/(log(x)^2+sigma^2);
fplot(P(x,0.1),[0.01 3],'DisplayName','P(x,0.1)','Color',[0 0.447 0.741],...
     'LineWidth',2);
hold on
fplot( P(x,0.5),[0.01 3],'DisplayName','P(x,0.5)','Color',[0.85 0.325 0.098],...
     'LineWidth',2);
fplot(P(x,1),[0.01 3],'DisplayName','P(x,1)','Color',[0.929 0.694 0.125],...
     'LineWidth',2);
axes1 = gca;
xlabel({'x'},'FontAngle','italic','FontName','Times New Roman');
title('The log-Cauchy distribution');
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);

% [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)
intvalue = simplify(int(P(xi,sigma),xi,[0,x]))

% [2 POINTS]
%
% Problem 4: Find the values of x that maximize P(x,0.5)

syms xatmax
xatmax = simplify(solve(diff(P(x,0.5),x)==0,x))

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

Patmax = subs(P(x,0.5),x,xatmax)
double(Patmax)   % Check which of the two solutions is greater
double(Patmax(2)) % The second one is the max...

% [2 POINTS]
%
%
%
%
% Problem 6: Solve the differential equation
% 6.1
% \frac{\textrm{d}N}{\textrm{dt}}=N_0 \left(1-N\left(\frac{c}{L}\right)\right)-\frac{N}{L}V_0 \left(1-N\;\frac{c}{L}\right)

clear all
syms N0 V0 c L t real
syms N(t)
eq = diff(N(t),t) == N0*(1-N(t)*c/L) - N(t)/L*V0*(1-N(t)*c/L)
IC = N(0)==0;
N(t) = simplify(dsolve(eq,IC,symvar(N(t))))

% [2 POINTS]
%
% 6.2 Plot N(t)

Nplot(t,V0) = subs(N(t),[L,c,N0],[500,3,1/10]);
figure
fplot(Nplot(t,11),[0 300],'DisplayName','V0=11','Color',[0 0.447 0.741],...
     'LineWidth',2);
hold on
fplot(Nplot(t,20),[0 300],'DisplayName','V0=20','Color',[0.85 0.325 0.098],...
     'LineWidth',2);
fplot(Nplot(t,45),[0 300],'DisplayName','V0=45','Color',[0.929 0.694 0.125],...
     'LineWidth',2);
axes1 = gca;
xlabel({'Time (sec)'},'FontAngle','italic','FontName','Times New Roman');
title('No. cars in road');
ylabel({'Number of cars'},'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]