I'm suppose to plot the conventional integral , but I'm getting wrong values

ok

[VBBV]

You need to specify the condition for the first 7 seconds given in your problem inside the loop as below

% Clear the workspace, close all figures, and clear command window
clc;
clear all;
close all;
% Given parameters
f1 = 3;     % Amplitude of the forcing function
f2 = 0;     % Second amplitude (set to 0)
T1 = 7;     % Time period for first interval
T2 = 30;    % Time period for second interval
T3 = 40;    % Time period for third interval
zeta = 0.1; % Damping ratio
wn = 3;     % Natural frequency
wd = wn * sqrt(1 - zeta^2); % Damped natural frequency
m = 1;      % Mass
% Define symbolic variables
syms t tau;
% Define the forcing function F1(t)
F1 = (f1*t/ T1);
% Calculate damping coefficient
c1 = (1 / (m * wd));
A = 1;
% Define the impulse response function x1(tau)
x1 = A*(exp(-zeta * (t - tau))) * (sin(wd * (t - tau)));
x2 = exp(-zeta * (t - tau)) * (sin(wd * (t - tau)));
% Define expressions for each interval of x(t)
x_1 = c1 * int(x1, tau, 0, t);
x_2 = c1 * (int(x1, tau, 0, T1) + int(0, tau, T1, t));
x_3 = c1 * (int(x1, tau, 0, T1) + int(0, tau, T1, T2) +  int(f2*x2, tau, T2, T3));
x_4 = c1 * (int(x1, tau, 0, T1) + int(0, tau, T1, T2) +  int(f2*x2, tau, T2, T3) + int(0, tau, T3, t));
% Total expression for x(t)
X_total = x_1 + x_2+ x_3 + x_4;
% Define time vector
t_values = linspace(0, T3, 1000);
% Evaluate x(t) for each time point
x_values = zeros(size(t_values));
for i = 1:length(t_values)
    if t_values(i) <= 7
        x_values(i) = double(subs(F1, t, t_values(i)));  % from the given conditions in your question 
    else
        x_values(i) = double(subs(X_total, t, t_values(i)));
    end
end
% Plot x(t)
plot(t_values, x_values, 'b', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('Displacement (m)');
title('Displacement Response x(t)');
grid on;

[Mohammad Adeeb]

thank you so much

Answer 1

[Torsten]

syms t x(t)
zeta = 0.1;
omegan = 3;
f1 = 1;
f2 = 0;
T1 = 7;
T2 = 30;
T3 = 40;
D2x = diff(x,t,2);
Dx = diff(x,t);
eqn = D2x + 2*zeta*omegan*Dx+omegan^2*x==f1*t/T1;
conds = [x(0)==0,Dx(0)==0];
x1 = dsolve(eqn,conds);
Dx1 = diff(x1,t);
eqn = D2x + 2*zeta*omegan*Dx+omegan^2*x==0;
conds = [x(T1)==subs(x1,t,T1),Dx(T1)==subs(Dx1,t,T1)];
x2 = dsolve(eqn,conds);
Dx2 = diff(x2,t);
eqn = D2x + 2*zeta*omegan*Dx+omegan^2*x==f2;
conds = [x(T2)==subs(x2,t,T2),Dx(T2)==subs(Dx2,t,T2)];
x3 = dsolve(eqn,conds);
Dx3 = diff(x3,t);
eqn = D2x + 2*zeta*omegan*Dx+omegan^2*x==0;
conds = [x(T3)==subs(x3,t,T3),Dx(T3)==subs(Dx3,t,T3)];
x4 = dsolve(eqn,conds);
x = piecewise(t>= 0 & t < T1,x1,t >=T1 & t < T2, x2, t>=T2 & t < T3,x3,t>=T3,x4);
fplot(x,[0 50])

[Mohammad Adeeb]

thank you so much for providing a help , does the first part seem to be correct ?

[Torsten]

What do you mean by "the first part" ?

[Mohammad Adeeb]

from zero to 7s approximatly

[Torsten]

You mean if the solution from 0 to 7 s in my code is correct ?

If the equation you want to solve for the damped spring/mass system is

x''+2*zeta*omega*x'+omega^2*x = f

I think it's correct.

But let's check:

syms t x(t)
zeta = 0.1;
omegan = 3;
f1 = 1;
f2 = 0;
T1 = 7;
T2 = 30;
T3 = 40;
D2x = diff(x,t,2);
Dx = diff(x,t);
eqn = D2x + 2*zeta*omegan*Dx+omegan^2*x==f1*t/T1;
conds = [x(0)==0,Dx(0)==0];
x1 = dsolve(eqn,conds);
Dx1 = diff(x1,t);
eqn = D2x + 2*zeta*omegan*Dx+omegan^2*x==0;
conds = [x(T1)==subs(x1,t,T1),Dx(T1)==subs(Dx1,t,T1)];
x2 = dsolve(eqn,conds);
Dx2 = diff(x2,t);
eqn = D2x + 2*zeta*omegan*Dx+omegan^2*x==f2;
conds = [x(T2)==subs(x2,t,T2),Dx(T2)==subs(Dx2,t,T2)];
x3 = dsolve(eqn,conds);
Dx3 = diff(x3,t);
eqn = D2x + 2*zeta*omegan*Dx+omegan^2*x==0;
conds = [x(T3)==subs(x3,t,T3),Dx(T3)==subs(Dx3,t,T3)];
x4 = dsolve(eqn,conds);
x = piecewise(t>= 0 & t < T1,x1,t >=T1 & t < T2, x2, t>=T2 & t < T3,x3,t>=T3,x4);
fplot(x,[0 50])
fun = @(t,x)[x(2);-2*zeta*omegan*x(2)-omegan^2*x(1)+f1*t/T1];
[T,Y] = ode45(fun,[0 T1],[0 0]);
hold on
plot(T,Y(:,1))
grid on

[Mohammad Adeeb]

thank you so much