Second order coupled differential equation
Afficher commentaires plus anciens
I need help in solving the second order coupled differential equations as I'm unable to get the right solution after spending many days.
I want to solve it with two for loops (One for frequency and other for Resistance), but unfortunatly it doesnt work well. I would be very happy if someone can help me in this regard.
I am sharing my code to be more precise.
X0 = [0; 0; 0]; % initial conditions
r = linspace(1, 50, 100); % resistance range
freq_vector = linspace(1, 20, 20); % frequency range
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
% Initialize array to store average power values over frequency range
avgpower_total = zeros(length(freq_vector), length(r));
% Loop over frequencies
for j = 1:length(freq_vector)
w = freq_vector(j); % Current frequency
avgpower = zeros(1, length(r));
tspan = linspace(0,100, 100);
F = 0.8*9.81*cos(w*tspan);
% Loop over resistance values
for i = 1:length(r)
R = r(i); % Current resistance
% Solve ODE using ode45
[t, X] = ode45(@(t, X) VoltResistFun2(t, X, R,F, tspan,w), tspan, X0, options);
% Extraction of voltage value from the solution
Voltage = X(:, 3); % Third state variable
% Calculation of power
P = Voltage.^2 / R;
% Integrated power over time
IntegratedPower = trapz(t, P);
% Calculate average power
avgpower(i) = IntegratedPower / (t(end) - t(1)); % Divide by total time
end
% Store average power for current frequency
avgpower_total(j, :) = avgpower;
% keyboard
end
% Integrate the power over the frequency range and divide by the range
avg_power_over_range = trapz(freq_vector, avgpower_total, 1) / (freq_vector(end) - freq_vector(1));
function dxdt = VoltResistFun2(t, X, R,F, tspan,w)
% Parameters
g = 9.81; % m/s^2
M = 0.01; % proof mass
t_p = 0.01; % thickness
A_p = 0.0001; % area
M_p = 0.00075; % proof mass
E_33 = 1.137e-8; % permittivity
k_33 = 0.75; % electromechanical coupling
e_33 = -1; % value of e_33
m = M + (1/3) * M_p;
C_p = E_33 * A_p / t_p;
theta = -(e_33 * A_p / t_p);
alpha = 1 ./ R;
w0 = 2 * pi * 20;
% F = 0.8 * g*cos(t); % external force
z = 0.02; % damping coefficient
f = interp1(tspan, F.*cos(w*t), t);
% Define state-space matrices
A = [0 1 0; -w0^2 -2*z*w0 -theta/m; 0 theta/C_p -alpha/C_p];
B = [0; -1; 0];
% Calculate derivative of state vector
% keyboard
dxdt = A\(X - B*f);
% Ensure dxdt is a column vector
% dxdt = dxdt(:);
end
Note :
eqution1 = x" +2zetaw_0x' +w_0^2x+theta/m= -F
F is defined
1/R = alpha
Thanks in advance
Réponse acceptée
X0 = [0; 0; 0]; % initial conditions
r = linspace(1, 50, 100); % resistance range
tspan = linspace(0,100, 100);
freq_vector = linspace(1, 20, 20); % frequency range
%options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
% Initialize array to store average power values over frequency range
avgpower_total = zeros(length(freq_vector), length(r));
% Loop over frequencies
for j = 1:length(freq_vector)
w = freq_vector(j); % Current frequency
avgpower = zeros(1, length(r));
F = @(t)0.8*9.81*cos(w*t);
% Loop over resistance values
for i = 1:length(r)
R = r(i); % Current resistance
% Solve ODE using ode45
[t, X] = ode45(@(t, X) VoltResistFun2(t, X, R,F), tspan, X0);%, options);
% Extraction of voltage value from the solution
Voltage = X(:, 3); % Third state variable
% Calculation of power
P = Voltage.^2 / R;
% Integrated power over time
IntegratedPower = trapz(t, P);
% Calculate average power
avgpower(i) = IntegratedPower / (t(end) - t(1)); % Divide by total time
end
% Store average power for current frequency
avgpower_total(j, :) = avgpower;
% keyboard
end
% Integrate the power over the frequency range and divide by the range
avg_power_over_range = trapz(freq_vector, avgpower_total, 1) / (freq_vector(end) - freq_vector(1))
plot(r,avg_power_over_range)
function dxdt = VoltResistFun2(t, X, R,F)
% Parameters
g = 9.81; % m/s^2
M = 0.01; % proof mass
t_p = 0.01; % thickness
A_p = 0.0001; % area
M_p = 0.00075; % proof mass
E_33 = 1.137e-8; % permittivity
k_33 = 0.75; % electromechanical coupling
e_33 = -1; % value of e_33
m = M + (1/3) * M_p;
C_p = E_33 * A_p / t_p;
theta = -(e_33 * A_p / t_p);
alpha = 1 ./ R;
w0 = 2 * pi * 20;
% F = 0.8 * g*cos(t); % external force
z = 0.02; % damping coefficient
f = F(t);
% Define state-space matrices
A = [0 1 0; -w0^2 -2*z*w0 -theta/m; 0 theta/C_p -alpha/C_p];
B = [0; -1; 0];
% Calculate derivative of state vector
% keyboard
dxdt = A\(X - B*f);
% Ensure dxdt is a column vector
% dxdt = dxdt(:);
end
8 commentaires
Muhammad
le 30 Mai 2024
Then you should check equations and parameters used. The code works technically - and that's where the forum could have helped.
According to your mathematical description, dxdt should be computed as
M = [1 0 0;0 m 0;0 0 Cp];
A = [0 1 0;-k -c theta;0 -theta 1/R_eq];
b = [0; -m*xb_doubledot;0];
dxdt = M\(A*X+b)
not as
dxdt = A\(X - B*f);
Steven Lord
le 30 Mai 2024
Rather than performing the division yourself in the ODE function, consider specifying M as the mass matrix using odeset. There are several examples that use a 'Mass' matrix listed in the Summary of ODE Examples and Files section of this documentation page. I can't tell from just the code whether or not any of those are closely related to the problem you're trying to solve, but if they are you could use that as a guide. [Even if they're not closely related, reading through one or two of them may still be useful.]
Muhammad
le 30 Mai 2024
Muhammad
le 30 Mai 2024
Hi @Muhammad
Your systems are technically stable by definition. However, their dominant eigenvalues have relatively large imaginary components. If you wish to obtain the desired response, then you must inject an external tunable signal that is derived from the desired response by applying some standard textbook formulas. I believe MATLAB also has the capability to automatically tune the tunable signal to produce the desired response.
R = 1:2:49;
for i = 1:numel(R)
% Parameters
g = 9.81; % m/s^2
M = 0.01; % proof mass
t_p = 0.01; % thickness
A_p = 0.0001; % area
M_p = 0.00075; % proof mass
E_33 = 1.137e-8; % permittivity
k_33 = 0.75; % electromechanical coupling
e_33 = -1; % value of e_33
m = M + (1/3) * M_p;
C_p = E_33 * A_p / t_p;
theta = -(e_33 * A_p / t_p);
alpha = 1 ./ R(i);
w0 = 2 * pi * 20;
% F = 0.8 * g*cos(t); % external force
z = 0.02; % damping coefficient
% f = F(t);
% Define state-space matrices
A = [0 1 0; -w0^2 -2*z*w0 -theta/m; 0 theta/C_p -alpha/C_p];
eigenvalues = eig(A)
end
%% Dominant Poles
format long
eigenvalues(1)
Muhammad
le 31 Mai 2024
Plus de réponses (0)
Catégories
En savoir plus sur Programming dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
