clear all
clc
% Define parameters
n = 3; % Number of degrees of freedom
m = eye(n); % Mass matrix (identity matrix for simplicity)
c = 0.1 * eye(n); % Damping matrix (proportional damping)
k = diag([11, 8, 9]); % Nonlinear stiffness matrix
% Time parameters
dt = 0.001; % Time step
t_final = 10; % Final time
t = 0:dt:t_final;
% Sine loading parameters
omega = 2*pi*0.5; % Frequency of the sine loading
F_sine = @(t) 200*10^3*sin(omega*t); % Sine loading function
% Initial conditions
x0 = zeros(n, 1); % Initial displacement
v0 = zeros(n, 1); % Initial velocity
% Newmark-Beta parameters
beta = 0.25; % Newmark-Beta parameter
gamma = 0.5; % Newmark-Beta parameter
% Initialize arrays
x = zeros(n, length(t)); % Displacement array
v = zeros(n, length(t)); % Velocity array
a = zeros(n, length(t)); % Acceleration array
% Initial conditions
x(:,1) = x0;
v(:,1) = v0;
a(:,1) = m \ (-c*v(:,1) - NonlinearStiffness(x(:,1), k));
% Newmark-Beta iteration
for i = 2:length(t)
% Predictor step
x_pred = x(:,i-1) + dt*v(:,i-1) + (dt^2/2) * ((1-2*beta)*a(:,i-1) + 2*beta*a(:,i-1));
v_pred = v(:,i-1) + dt*((1-gamma)*a(:,i-1));
% Corrector step
a_pred = m \ (-c*v_pred - NonlinearStiffness(x_pred, k) + F_sine(t(i)));
x(:,i) = x(:,i-1) + dt*v_pred + (dt^2/2) * (beta*a_pred + (1-beta)*a(:,i-1));
v(:,i) = v_pred + dt*((1-gamma)*a_pred + dt*gamma*a(:,i));
% Update acceleration for next iteration
a(:,i) = a_pred;
end
% Plot results
plot(t, x(1:3,:));
xlabel('Time');
ylabel('Displacement DOF 1 ')
xlim([0 5]); grid
% Nonlinear stiffness function
function F = NonlinearStiffness(x, k)
F = k * x.^3; % Element-wise cubic nonlinear stiffness
end
