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
