Michael Chuah's page for MAS.864 The Nature of Mathematical Modelling

close all; clear all; clc;

%% 3-1(c)
m = 1; k = 1; b = 0.1;
syms w;

magC = ((k-m.*w.^2).^2+(w.*b).^2).^(-1/2);
phaseC = atan(w.*b./(k-m.*w.^2));

freq = -pi:0.01:pi;
figure; 
subplot(2,1,1); plot(freq,subs(magC,freq)); grid on; 
title('Magnitude plot'); xlabel('Frequency (rad)'); ylabel('Magnitude'); ylim([0 10]);
subplot(2,1,2); plot(freq,subs(phaseC,freq)); grid on;
title('Phase plot'); xlabel('Frequency (rad)'); ylabel('Phase (rad)'); ylim([-pi/2 pi/2]);
saveas(gcf,'3-1(c) Magnitude and phase plot.png');
saveas(gcf,'3-1(c) Magnitude and phase plot.fig');

%% 3-1(d)
clc;
syms m k b w;
S = solve(m^3*w^4+(m*b^2-2*k*m^2)*w^2+m*k^2-2*k*b^2,w)
pretty(S)

w1 =  (-(b*(4*k*m)^(1/2) - 2*k*m)/(2*m^2))^(1/2);
w2 = ((2*k*m + b*(4*k*m)^(1/2))/(2*m^2))^(1/2);

pretty(w1)
pretty(w2)

pretty(subs(diff(w1,b),b,0))

pretty(diff((1-exp(-b/((m*k)^(1/2))))^(-1),b))

S =
 
  (-(b*(b^2 + 4*k*m)^(1/2) - 2*k*m + b^2)/(2*m^2))^(1/2)
   ((2*k*m + b*(b^2 + 4*k*m)^(1/2) - b^2)/(2*m^2))^(1/2)
 -(-(b*(b^2 + 4*k*m)^(1/2) - 2*k*m + b^2)/(2*m^2))^(1/2)
  -((2*k*m + b*(b^2 + 4*k*m)^(1/2) - b^2)/(2*m^2))^(1/2)
 
 
  +-     -+ 
  |   #1  | 
  |       | 
  |   #2  | 
  |       | 
  |  -#1  | 
  |       | 
  |  -#2  | 
  +-     -+ 
   
  where 
   
           /       2         1/2            2 \1/2 
           |   b (b  + 4 k m)    - 2 k m + b  | 
     #1 == | - ------------------------------ | 
           |                   2              | 
           \                2 m               / 
   
           /             2         1/2    2 \1/2 
           | 2 k m + b (b  + 4 k m)    - b  | 
     #2 == | ------------------------------ | 
           |                 2              | 
           \              2 m               /
 
  /                  1/2 \1/2 
  | 2 k m - b (4 k m)    | 
  | -------------------- | 
  |            2         | 
  \         2 m          /
 
  /                  1/2 \1/2 
  | 2 k m + b (4 k m)    | 
  | -------------------- | 
  |            2         | 
  \         2 m          /
 
            1/2 
     (4 k m) 
  - ------------- 
       2 / k \1/2 
    4 m  | - | 
         \ m /
 
                /      b     \ 
             exp| - -------- | 
                |        1/2 | 
                \   (k m)    / 
  - ----------------------------------- 
    /    /      b     \     \2      1/2 
    | exp| - -------- | - 1 |  (k m) 
    |    |        1/2 |     | 
    \    \   (k m)    /     /
>>