I'm suppose to plot the response of the system but i'm getting wrong values

i was trying to solve the system , but i'm getting wrong va;ues , what am i missing?here is my code

clc;
clear all;
close all;
% Define system parameters
params = struct();
% System parameters
params.M1 = 6e-3;      % tip Mass of beam I in kg
params.M2 = 5e-3;      % tip Mass of beam II in kg
params.L1 = 7e-2;      % Length of beam I in meters
params.L2 = 4e-2;      % Length of beam II in meters
params.w1 = 2e-2;      % Width of beam I in meters
params.w2 = 2e-2;      % Width of beam II in meters
params.h1 = 1e-3;      % Thickness of beam I in meters
params.h2 = 1e-3;      % Thickness of beam II in meters
params.E1 = 69e9;      % Young's modulus of beam I in Pascals
params.E2 = 2.344e9;   % Young's modulus of beam II in Pascals
params.c1 = 0.1;       % Damping of beam I in Ns/m
params.c2 = 0.1;       % Damping of beam II in Ns/m
params.c12 = 0;        % Couplind damping in Ns/m
params.c21 = 0;        % Couplind damping in Ns/m
params.theta = 0.4e-3; % Piezoelectric constant in m/V
params.R = 10e6;       % Resistance in ohms
params.Cp = 72e-9;     % Piezoelectric capacity in farads
params.F0 = 5;         % Amplitude of base excitation force in Newton
params.omega_a = 10;   % Angular frequency of base excitation in rad/s
params.epsilon_0 = 8.854e-12; % Vacuum permittivity in F/m
params.tau = 1e-3;     % Time constant in seconds
params.epsilon_r = 1;  % Relative permittivity
params.do_2 = 1e-3;    % Constant
params.sigma = 1e-9;   % Conductivity in S/m
params.rho1 = 2700;    % Density for Beam I in kg/m^3
params.rho2 = 1220;    % Density for Beam II in kg/m^3
params.lt = 1e-2;      % the length of the Tribo electric material in m
params.wt = 1e-2;      % the width of the Tribo electric material in m
params.S = params.lt * params.wt; % contact area of the tribo material
params.gi = 6e-1;      % the gap in m
params.I1 = (1/12) * params.w1 * (params.h1^3); % moment of inertia for first beam 
params.I2 = (1/12) * params.w2 * (params.h2^3); % moment of inertia for second beam
params.alpha = 2*params.epsilon_0*params.epsilon_r*params.S;
params.alpha2 = params.epsilon_0*params.R*params.S;
params.pt = 2e-2 ; %thickness of pizo
params.sigma = 1e-6; %surfce charge density for PDMS
params.gap = 6e-3; % the gap For PDMS 
params.alpha3 = ((params.pt/params.epsilon_0)+params.do_2-params.gap)+((params.sigma/(params.R*params.epsilon_0))*(params.do_2-params.gap));
% Calculate equivalent masses
params.m1 = params.rho1 * params.L1 * params.w1 * params.h1; 
params.m2 = params.rho2 * params.L2 * params.w2 * params.h2;
params.Meq1 = params.M1 + (33 * params.m1 * params.L1) / 140;
params.Meq2 = params.M2 + (33 * params.m2 * params.L2) / 140;
% Calculate ko
ko = 6 * params.E1 * params.I1 / ((4 * params.E1 * params.I1 * params.L1^3 * params.L2^3) + (3 * params.E2 * params.I2 * params.L1^4 * params.L2^2));
% Calculate K matrix
K = ko * [2 * ((2 * params.E1 * params.I1 * params.L2^3) + 6 * params.E2 * params.I2 * (params.L1 * params.L2^2 - params.L1^2 * params.L2) + (2 * params.E2 * params.I2 * params.L1^3)), ...
    (3 * params.E2 * params.I2 * params.L1^2 * params.L2) - (2 * params.E2 * params.I2 * params.L1^3);
    (3 * params.E2 * params.I2 * params.L1^2 * params.L2) - (2 * params.E2 * params.I2 * params.L1^3), ...
    2 * params.E2 * params.I2 * params.L1^3];
% Extract individual elements
params.k11 = K(1, 1);
params.k12 = K(1, 2);
params.k21 = K(2, 1);
params.k22 = K(2, 2);
% Define state-space matrices
A = [0, 1, 0, 0, 0, 0;
    -(params.k11/params.Meq1), -(params.c1/params.Meq1), -(params.k12/params.Meq1), 0, 0, -(params.theta/params.Meq1);
    0, 0, 0, 1, 0, 0;
    -(params.k21/params.Meq2), 0, -(params.k22/params.Meq2), -(params.c21/params.Meq2), 0, 0;
    0, 0, 0, 0, 0, 1/params.R;
    0, 0, 0, 0, 0, -(params.theta/params.Cp)];
% Define time span
tspan = [0 20]; % Time span from 0 to 10 seconds
% Define initial conditions
x0 = [1;1;1;1;1;1];
% Define input function (if any)
a = @(t) params.F0*sin(params.omega_a*t); % No input function for now
% Define the function for the system of first-order differential equations
dxdt = @(t, x) [x(2);
    -(params.c1/params.Meq1)*x(2) - (params.c12/params.Meq1)*x(4) - (params.k11/params.Meq1)*x(1) - (params.k12/params.Meq1)*x(3) - (params.theta/params.Meq1)*x(6) - a(t);
    x(4);
    -(params.c2/params.Meq2)*x(4) - (params.c21/params.Meq2)*x(2) - (params.k22/params.Meq2)*x(3) - (params.k21/params.Meq2)*x(1) - ((x(5))^2)/(params.alpha)- a(t);
    ((x(5))/params.alpha2) * (params.alpha3 + params.do_2 - params.gap) + params.alpha3 * (params.do_2 - params.gap);
    x(6)/params.R + params.Cp * x(6)];
% Solve the system using ode45
[t, x] = ode45(dxdt, tspan, x0);
% Extracting the results
y1 = x(:,1);
y1dot = x(:,2);
y2 = x(:,3);
y2dot = x(:,4);
q = x(:,5);
v = x(:,6);
% Plot the results if needed
% For example, plot y1 and y2 over time
plot(t, q, 'b', t, v, 'r');
xlabel('Time');
ylabel('Displacement');
legend('y1', 'y2');
title('Displacement Response');

This the equations :

[Torsten]

Why does params.gap equal x(3) in equation (5) ?

And note that in equation (5), you divide by params.alpha2 which is in the order of 1e-12. This will blow up the time derivative tremendously.

And you will have to define your dxdt in a function instead of a function handle because the differential equation for y2 depends on a condition.

[Mohammad Adeeb]

alpha is a constant , and this is the value of it , params.gap = 6e-6 not x3

and for dxdt , what do you mean by a function , you mean seperate one ? if so , can you please guide me through the process of definning it and call to the ode45

Answer 1

[Torsten]

alpha is a constant , and this is the value of it ,

I doubt that the rate of change at t=0 for x(5) should be 2.55e17 as you get when you call dxdt for the vector of initial values (see below). The main contribution is a factor of 1e12 caused by params.alpha2.

params.gap = 6e-6 not x3

The equation for q2 in your mathematical description is written with y2, not with params.gap.

and for dxdt , what do you mean by a function , you mean seperate one ? if so , can you please guide me through the process of definning it and call to the ode45

By defining dxdt in a function file as below, you can easily use an if-condition to define the equation for y2 state-dependent:

% Define time span
tspan = [0 20]; % Time span from 0 to 10 seconds
% Define initial conditions
x0 = [1;1;1;1;1;1];
% Solve the system using ode45
fun(0,x0)

ans = 6×1

1.0e+17 * 0.0000 -0.0000 0.0000 -0.0056 2.5512 0.0000

[t, x] = ode45(@fun, tspan, x0);
% Extracting the results
y1 = x(:,1);
y1dot = x(:,2);
y2 = x(:,3);
y2dot = x(:,4);
q = x(:,5);
v = x(:,6);
% Plot the results if needed
% For example, plot y1 and y2 over time
plot(t, q, 'b', t, v, 'r');
xlabel('Time');
ylabel('Displacement');
legend('y1', 'y2');
title('Displacement Response');

function dxdt = fun(t,x)
% Define system parameters
params = struct();
% System parameters
params.M1 = 6e-3;      % tip Mass of beam I in kg
params.M2 = 5e-3;      % tip Mass of beam II in kg
params.L1 = 7e-2;      % Length of beam I in meters
params.L2 = 4e-2;      % Length of beam II in meters
params.w1 = 2e-2;      % Width of beam I in meters
params.w2 = 2e-2;      % Width of beam II in meters
params.h1 = 1e-3;      % Thickness of beam I in meters
params.h2 = 1e-3;      % Thickness of beam II in meters
params.E1 = 69e9;      % Young's modulus of beam I in Pascals
params.E2 = 2.344e9;   % Young's modulus of beam II in Pascals
params.c1 = 0.1;       % Damping of beam I in Ns/m
params.c2 = 0.1;       % Damping of beam II in Ns/m
params.c12 = 0;        % Couplind damping in Ns/m
params.c21 = 0;        % Couplind damping in Ns/m
params.theta = 0.4e-3; % Piezoelectric constant in m/V
params.R = 10e6;       % Resistance in ohms
params.Cp = 72e-9;     % Piezoelectric capacity in farads
params.F0 = 5;         % Amplitude of base excitation force in Newton
params.omega_a = 10;   % Angular frequency of base excitation in rad/s
params.epsilon_0 = 8.854e-12; % Vacuum permittivity in F/m
params.tau = 1e-3;     % Time constant in seconds
params.epsilon_r = 1;  % Relative permittivity
params.do_2 = 1e-3;    % Constant
params.sigma = 1e-9;   % Conductivity in S/m
params.rho1 = 2700;    % Density for Beam I in kg/m^3
params.rho2 = 1220;    % Density for Beam II in kg/m^3
params.lt = 1e-2;      % the length of the Tribo electric material in m
params.wt = 1e-2;      % the width of the Tribo electric material in m
params.S = params.lt * params.wt; % contact area of the tribo material
params.gi = 6e-1;      % the gap in m
params.I1 = (1/12) * params.w1 * (params.h1^3); % moment of inertia for first beam 
params.I2 = (1/12) * params.w2 * (params.h2^3); % moment of inertia for second beam
params.alpha = 2*params.epsilon_0*params.epsilon_r*params.S;
params.alpha2 = params.epsilon_0*params.R*params.S;
params.pt = 2e-2 ; %thickness of pizo
params.sigma = 1e-6; %surfce charge density for PDMS
params.gap = 6e-3; % the gap For PDMS 
params.alpha3 = ((params.pt/params.epsilon_0)+params.do_2-params.gap)+((params.sigma/(params.R*params.epsilon_0))*(params.do_2-params.gap));
% Calculate equivalent masses
params.m1 = params.rho1 * params.L1 * params.w1 * params.h1; 
params.m2 = params.rho2 * params.L2 * params.w2 * params.h2;
params.Meq1 = params.M1 + (33 * params.m1 * params.L1) / 140;
params.Meq2 = params.M2 + (33 * params.m2 * params.L2) / 140;
% Calculate ko
ko = 6 * params.E1 * params.I1 / ((4 * params.E1 * params.I1 * params.L1^3 * params.L2^3) + (3 * params.E2 * params.I2 * params.L1^4 * params.L2^2));
% Calculate K matrix
K = ko * [2 * ((2 * params.E1 * params.I1 * params.L2^3) + 6 * params.E2 * params.I2 * (params.L1 * params.L2^2 - params.L1^2 * params.L2) + (2 * params.E2 * params.I2 * params.L1^3)), ...
    (3 * params.E2 * params.I2 * params.L1^2 * params.L2) - (2 * params.E2 * params.I2 * params.L1^3);
    (3 * params.E2 * params.I2 * params.L1^2 * params.L2) - (2 * params.E2 * params.I2 * params.L1^3), ...
    2 * params.E2 * params.I2 * params.L1^3];
% Extract individual elements
params.k11 = K(1, 1);
params.k12 = K(1, 2);
params.k21 = K(2, 1);
params.k22 = K(2, 2);
% Define state-space matrices
A = [0, 1, 0, 0, 0, 0;
    -(params.k11/params.Meq1), -(params.c1/params.Meq1), -(params.k12/params.Meq1), 0, 0, -(params.theta/params.Meq1);
    0, 0, 0, 1, 0, 0;
    -(params.k21/params.Meq2), 0, -(params.k22/params.Meq2), -(params.c21/params.Meq2), 0, 0;
    0, 0, 0, 0, 0, 1/params.R;
    0, 0, 0, 0, 0, -(params.theta/params.Cp)];
% Define input function (if any)
a =  params.F0*sin(params.omega_a*t); % No input function for now
% Define the function for the system of first-order differential equations
dxdt = [x(2);
    -(params.c1/params.Meq1)*x(2) - (params.c12/params.Meq1)*x(4) - (params.k11/params.Meq1)*x(1) - (params.k12/params.Meq1)*x(3) - (params.theta/params.Meq1)*x(6) - a;
    x(4);
    -(params.c2/params.Meq2)*x(4) - (params.c21/params.Meq2)*x(2) - (params.k22/params.Meq2)*x(3) - (params.k21/params.Meq2)*x(1) - ((x(5))^2)/(params.alpha)- a;
    ((x(5))/params.alpha2) * (params.alpha3 + params.do_2 - params.gap) + params.alpha3 * (params.do_2 - params.gap);
    x(6)/params.R + params.Cp * x(6)];
end

[Mohammad Adeeb]

ok , i got your point , but the thing is I'm still getting wrong output , the output should have 2 peaks , i for y1 and the other is for y2 , i don't know I'm doing things wrong , the peaks is for omega 1and2 , plus the q2 dot should have a value when plotting it with time

thanks alot for the help

[Torsten]

We cannot help with your model and the output it gives. We can only show you potential sources of errors in your code.

[Sam Chak]

@Mohammad Adeeb, I noticed that the state matrix in your code is unused. Could you please clarify this? Additionally, have you verified if the differential equations are correctly written in the code? If possible, it would be helpful if you could provide an image of the expected graph for and .

[Sam Chak]

Hi @Mohammad Adeeb,

I have some concerns about the reliability of the differential equations presented in your image. Typically, reputable journals or publishers adhere to standard mathematical notations. It appears that the equations might have been typed using Microsoft Equation Editor or LaTeX without being thoroughly proofread by professionals.

If the differential equations in your image are indeed correct, they can be rearranged as follows: