Can anyone help me to get the graph of x2 vs t and x1 vs t using ode 45? md^2x2/dt^​2-k(x1-x2)​+c(dx2/dt)​/h=0, m=1,k=100,​c=0.001,h=​0.01, at t=0,x1=0,x​2=0,dx1/dt​=0.01 ,dx2/dt=0

[Torsten]

How is x1(t) defined ?

Answer 1

[Simar]

Hi Puja,

I understand that you are facing an issue in obtaining the graph of two functions, x2 versus t and x1 versus t, by solving a second-order ordinary differential equation (ODE) using MATLAB's “ode45” solver. Please refer to the solution provided below.

Firstly, create a separate function file since the system of ODEs provided is a bit complex, it is better to create a separate function file rather than an anonymous function. Create a new file named “odesys.m” as follow:

function dxdt = odesys(t, x, m, k, c, h) 
% Unpack the state vector 
x1 = x(1); 
dx1dt = x(2); 
x2 = x(3); 
dx2dt = x(4); 
% Define the system of first-order ODEs 
d2x2dt2 = k * (x1 - x2) / m - (c * dx2dt) / (m * h); 
dxdt = [dx1dt; 0; dx2dt; d2x2dt2]; % The derivative of x1 is dx1dt, and its second derivative is 0 since it is not defined in the problem 
end 

The MATLAB code provided defines and solves a system of ordinary differential equations (ODEs) representing a mechanical system with two degrees of freedom, x1 and x2. The system includes parameters for mass (m), spring constant (k), damping coefficient (c), and another damping-related constant (h). The “odesys” function specifies how the system evolves over time by defining the relationships between the displacements and velocities of x1 and x2.

The ODE solver “ode45” is used to numerically integrate the system of equations over a specified time span, starting from given initial conditions for the displacements and velocities of x1 and x2. The solver computes the values of x1 and x2 at various points in time, which can then be plotted to visualize the behaviour of the system.

Afterwards in the main script follow the steps:

1. Define the parameters that will be used in the differential equation. m is the mass, k is the spring constant, c is the damping coefficient, and h is a constant that appears in the damping term.

m = 1; 
k = 100; 
c = 0.001; 
h = 0.01; 

2. Set initial conditions for the system. x1_0 and x2_0 are the initial displacements of x1 and x2, respectively, and dx1dt_0 and dx2dt_0 are the initial velocities.

x1_0 = 0; 
x2_0 = 0; 
dx1dt_0 = 0.01; 
dx2dt_0 = 0; 

3. The initial_conditions array combines the initial conditions into a single column vector. It will be used as the starting point for the “ode45” solver.

initial_conditions = [x1_0; dx1dt_0; x2_0; dx2dt_0]; 

4. tspan is a vector specifying the start and end times for the simulation. In this case, it starts at t = 0 and ends at t = 10. Adjust the end time based on how long one wants to simulate the system.

tspan = [0 10]; 

5. The odesys function (defined in odesys.m) created in beginning. This function file “odesys.m” defines the system of first-order ordinary differential equations (ODEs) that ode45 will solve.

function dxdt = odesys(t, x, m, k, c, h) 
... 
    End 

6. Solving the ODE using ode45. The solver returns a vector t of time points and a matrix x where each row corresponds to the solution at a time point in t.

[t, x] = ode45(@(t, x) odesys(t, x, m, k, c, h), tspan, initial_conditions); 

7. Extracting x1 and x2 from the solution. The solution matrix x has four columns, corresponding to x1, dx1dt, x2, and dx2dt. Since interested in the displacements x1 and x2, so extracting the first and third columns.

x1 = x(:, 1); 
x2 = x(:, 3); 

8. Plotting the results and generating two separate figures: one for x2 vs t and another for x1 vs t. Each figure opens a new window, plots the respective displacement against time, and labels the axes and titles appropriately.

% Plot x2 vs t 
figure; 
plot(t, x2); 
title('x2 vs t'); 
xlabel('Time, t'); 
ylabel('Displacement, x2'); 
% Plot x1 vs t 
figure; 
plot(t, x1); 
title('x1 vs t'); 
xlabel('Time, t'); 
ylabel('Displacement, x1'); 

Here is the compiled code for main script:

% Define the parameters 
m = 1; 
k = 100; 
c = 0.001; 
h = 0.01; 
% Initial conditions 
x1_0 = 0; 
x2_0 = 0; 
dx1dt_0 = 0.01; 
dx2dt_0 = 0; 
% Initial state vector 
initial_conditions = [x1_0; dx1dt_0; x2_0; dx2dt_0]; 
% Time span for the simulation 
tspan = [0 10]; % You can adjust the end time as needed 
% Solve the ODE using ode45 
[t, x] = ode45(@(t, x) odesys(t, x, m, k, c, h), tspan, initial_conditions); 
% Extract x1 and x2 from the solution 
x1 = x(:, 1); 
x2 = x(:, 3); 
% Plot x2 vs t 
figure; 
plot(t, x2); 
title('x2 vs t'); 
xlabel('Time, t'); 
ylabel('Displacement, x2'); 
% Plot x1 vs t 
figure; 
plot(t, x1); 
title('x1 vs t'); 
xlabel('Time, t'); 
ylabel('Displacement, x1'); 

Following are the plots:

Hope it helps!

Best Regards,

Simar