Write a MATLAB script for plotting (a) the non-dimensional response magnitude for a system
with harmonically moving base shown in Fig. the response phase angle for system with
harmonically moving base
zeta= [0.05; 0.1; 0.15; 0.25; 0.5; 1.25; 1.5]; % damping factors
r= [0:0.01:3]; %frequency ratio
for k=1: length (zeta)
G(k,:)=sqrt((1+(2*zeta(k)*r).^2)./((1–r.^2).^2+(2*zeta(k)*r).^2));
phi(k,:)=atan2(2*zeta(k)*r.^3,1–r.^2+(2*zeta(k)*r).^2); end
figure (1)
plot(r, G)
xlabel (‘\omega/\omega_n’)
ylabel (‘|x (i\omega)|/A’)
grid
legend(‘\zeta_1=0.05’,‘\zeta_2=0.1’,‘\zeta_3=0.15’,‘\zeta_4=0.25’,‘\zeta_5=
0.5’, ‘\zeta_6=1.25’,‘\zeta_7=1.5’)
figure (2)
plot(r, phi)
xlabel (‘\omega/\omega_n’)
ylabel (‘\phi (\omega)’)
grid
ha=gca;
set (ha,’ytick’,[0:pi/2:pi])
set(ha,’yticklabel’,{[];’pi/2';‘p’})
legend(‘\zeta_1=0.05’,‘\zeta_2=0.1’,‘\zeta_3=0.15’,‘\zeta_4=0.25’,‘\zeta_5=
0.5’, ‘\zeta_6=1.25’, ‘\zeta_7=1.5’)
% Numerical Integration and Differentiation
% Implements Trapezoidal Rule, Simpson's Rule, and Finite Difference Differentiation
% Sample function: f(x) = x^2 + 2x + 1
f = @(x) x.^2 + 2*x + 1;
% Trapezoidal Rule for Integration
function I = trapezoidal_rule(f, a, b, n)
% f: function handle
% a, b: integration limits
% n: number of subintervals
h = (b - a) / n;
x = a:h:b;
y = f(x);
I = (h/2) * (y(1) + 2*sum(y(2:end-1)) + y(end));
end
% Simpson's Rule for Integration
function I = simpson_rule(f, a, b, n)
% f: function handle
% a, b: integration limits
% n: number of subintervals (must be even)
if mod(n, 2) ~= 0
error('Number of subintervals must be even for Simpson''s Rule');
end
h = (b - a) / n;
x = a:h:b;
y = f(x);
I = (h/3) * (y(1) + 4*sum(y(2:2:end-1)) + 2*sum(y(3:2:end-2)) +
y(end));
End
% Finite Difference for Differentiation (Central Difference)
function df = central_difference(f, x, h)
% f: function handle
% x: point at which to compute derivative
% h: step size
df = (f(x + h) - f(x - h)) / (2*h);
end
% Main script to demonstrate the methods
% Integration parameters
a = 0; % Lower limit
b = 2; % Upper limit
n = 100; % Number of subintervals (even for Simpson's Rule)
% Compute integrals
trap_result = trapezoidal_rule(f, a, b, n);
simp_result = simpson_rule(f, a, b, n);
% Differentiation parameters
x_point = 1; % Point to compute derivative
h = 0.01; % Step size for differentiation
% Compute derivative
diff_result = central_difference(f, x_point, h);
% Display results
fprintf('Trapezoidal Rule Integral (from %d to %d): %.6f\n', a, b,
trap_result);
fprintf('Simpson''s Rule Integral (from %d to %d): %.6f\n', a, b,
simp_result);
fprintf('Derivative at x = %.2f: %.6f\n', x_point, diff_result);
% Compare with analytical results
% Analytical integral of x^2 + 2x + 1 from 0 to 2: [x^3/3 + x^2 + x] = 14/3
% Analytical derivative at x = 1: 2x + 2 = 4
analytical_integral = 14/3;
analytical_derivative = 4;
fprintf('Analytical Integral: %.6f\n', analytical_integral);
fprintf('Analytical Derivative at x = %.2f: %.6f\n', x_point,
analytical_derivative);
% Plot the function
x = linspace(a, b, 100);
y = f(x);
figure;
plot(x, y, 'b-', 'LineWidth', 2);
title('Function f(x) = x^2 + 2x + 1');
xlabel('x');
ylabel('f(x)');
grid on;
% MATLAB Code for Engineering Analysis
% Covers Basic Engineering Mechanics, Mechanical Vibrations, Control Systems,
% Statistics, and Circuit Dynamics
%% 1. Basic Engineering Mechanics: Beam Deflection
% Simply supported beam with a point load at the center
L = 2; % Beam length (m)
P = 1000; % Point load (N)
E = 200e9; % Young's modulus (Pa)
I = 1e-6; % Moment of inertia (m^4)
x = linspace(0, L, 100); % Position along beam
deflection = @(x) (-P * x .* (L^2 - x.^2)) / (48 * E * I);
% Deflection formula
y = deflection(x);
figure(1);
plot(x, y * 1000, 'b-', 'LineWidth', 2);
% Convert to mm for visualization
title('Beam Deflection under Point Load');
xlabel('Position along beam (m)');
ylabel('Deflection (mm)');
grid on;
%% 2. Mechanical Vibrations: Single DOF System
% Free vibration of a mass-spring damper system
m = 1; % Mass (kg)
k = 100; % Spring constant (N/m)
c = 2; % Damping coefficient (Ns/m)
omega_n = sqrt(k/m); % Natural frequency
zeta = c / (2 * sqrt(m * k)); % Damping ratio
t = 0:0.01:5; % Time vector
x0 = 0.1; % Initial displacement (m)
v0 = 0; % Initial velocity (m/s)
% Underdamped response
x = exp(-zeta * omega_n * t).* (x0 * cos(omega_n * sqrt(1 - zeta^2) * t) +
...
(v0 + zeta * omega_n * x0) / (omega_n * sqrt(1 - zeta^2)) * sin(omega_n
* sqrt(1 - zeta^2) * t));
figure(2);
plot(t, x, 'r-', 'LineWidth', 2);
title('Free Vibration of Mass-Spring-Damper System');
xlabel('Time (s)');
ylabel('Displacement (m)');
grid on;
%% 3. Control Systems: Step Response of a Second-Order System
% Transfer function: G(s) = omega_n^2 / (s^2 + 2*zeta*omega_n*s +
omega_n^2)
num = omega_n^2; % Numerator
den = [1, 2 * zeta * omega_n, omega_n^2]; % Denominator
sys = tf(num, den); % Transfer function
t = 0:0.01:5;
[y, t] = step(sys, t); % Step response
figure(3);
plot(t, y, 'g-', 'LineWidth', 2);
title('Step Response of Second-Order System');
xlabel('Time (s)');
ylabel('Output');
grid on;
%% 4. Statistics: Analysis of Random Data
% Generate random data for stress in a material (MPa)
rng(0); % Set seed for reproducibility
stress = 100 + 20 * randn(1000, 1); % Normal distribution
mean_stress = mean(stress);
std_stress = std(stress);
% Histogram
figure(4);
histogram(stress, 20, 'FaceColor', 'b');
title('Stress Distribution in Material');
xlabel('Stress (MPa)');
ylabel('Frequency');
grid on;
fprintf('Mean Stress: %.2f MPa\n', mean_stress);
fprintf('Standard Deviation: %.2f MPa\n', std_stress);
%% 5. Dynamics of Electrical Circuits: RLC Circuit Response
% Series RLC circuit with step input
R = 10; % Resistance (ohms)
L = 0.1; % Inductance (H)
C = 1e-3; % Capacitance (F)
num = [1/(L*C)]; % Transfer function numerator (Vout/Vin)
den = [1, R/L, 1/(L*C)]; % Denominator
sys_rlc = tf(num, den);
t = 0:0.001:0.5;
[y_rlc, t] = step(sys_rlc, t); % Step response
figure(5);
plot(t, y_rlc, 'm-', 'LineWidth', 2);
title('RLC Circuit Step Response');
xlabel('Time (s)');
ylabel('Output Voltage (V)');
grid on;
% Display key results
fprintf('Beam Max Deflection: %.2f mm\n', max(abs(y)) * 1000);
fprintf('Vibration Natural Frequency: %.2f rad/s\n', omega_n);
fprintf('Damping Ratio: %.2f\n', zeta);