Second order coupled differential equation
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 advanceI 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 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 ode45, second order coupled equations, energy harvester, avg power calculation MATLAB Answers — New Questions
​
