MATLAB: Empty sym 0-by-1

empty sym

Hello,
I followed a few responses to the empty sym 0-by-1 issue and they pointed out that I had to indicate which parameter I was solving for.
I had already specified that in my code. I am trying to solve for R_L for the maximum of the plotted function, the "R_l_optimal" function is what is giving me grief. I know it looks around 70, but I want to calculate that specifically, and I keep getting the empty sym 0-by-1 result.
Can you see where my issue is?
Thank you!
clear all
syms omega t R_l s A_in Y
%Inputs
A_o = 3;
R_c = 40;
L_c = 0.051;
f=186;
omega_in = 2*pi*f;
A = A_o*sin(omega*t);
M=0.01;
%Theta from part a)
theta = 0.244;
%Spring coefficient (N/m)
k = 13660;
%Damping coefficient (N*s/m)
C_d = 0.07;
%Impedance equations
Z1 = R_c+R_l+L_c*s;
Z2 = C_d+k/s+M*s;
%From part a) transfer function
%Output voltage V_l
V_l = (theta*(Y*s)*R_l)/Z1;
%Input excitation force Z_in
A_in=(Y*s*Z2+theta*V_l/R_l)/(M);
%output voltage ratio to input force
Transfer_func = simplify(V_l/A_in);
Transfer_func(s) = vpa(simplify(V_l/A_in),5);
Transfer_func(omega) = subs(Transfer_func, {s},{1j*omega});
amplitude_TF = abs(Transfer_func(omega_in));
angle_TF = angle(Transfer_func(omega_in));
V_l_steady_state = amplitude_TF*(sin(omega_in*t+angle_TF));
P_ave_output = (amplitude_TF/(sqrt(2)))^2/R_l;
P_ave_diff = diff(P_ave_output);
P_ave_max = P_ave_diff == 0;
R_l_optimal = (solve(P_ave_max, R_l))
fplot(P_ave_output,[10 200])
xticks([10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200])
xline(70)
ylabel('Average Power (Watts)')
xlabel('Load Resistance R_l (Ohms)')

Best Answer

The problem is one that you couldn’t have anticipated, and that it took me a while to discover. It apparently has to do with frailties in the solve function.
In addition to requiring ‘R_l’ to be , your code needs one extra line: 
P_ave_diff = simplify(P_ave_diff, 'Steps',500)                      % <— Add This


and then: 
R_l_optimal =
72.275099696632748765256742994398
The Full (Slightly Revised) Code — 
syms omega t R_l s A_in Y
assume(R_l >= 0)                                                    % <— Add This
%Inputs
A_o = 3;
R_c = 40;
L_c = 0.051;
f=186;
omega_in = 2*pi*f;
A = A_o*sin(omega*t);
M=0.01;
%Theta from part a)
theta = 0.244;
%Spring coefficient (N/m)
k = 13660;
%Damping coefficient (N*s/m)
C_d = 0.07;
%Impedance equations
Z1 = R_c+R_l+L_c*s;
Z2 = C_d+k/s+M*s;
%From part a) transfer function
%Output voltage V_l
V_l = (theta*(Y*s)*R_l)/Z1;
%Input excitation force Z_in
A_in=(Y*s*Z2+theta*V_l/R_l)/(M);
%output voltage ratio to input force
Transfer_func = simplify(V_l/A_in);
Transfer_func(s) = vpa(simplify(V_l/A_in),5);
Transfer_func(omega) = subs(Transfer_func, {s},{1j*omega});
amplitude_TF = abs(Transfer_func(omega_in));
angle_TF = angle(Transfer_func(omega_in));
V_l_steady_state = amplitude_TF*(sin(omega_in*t+angle_TF));
P_ave_output = (amplitude_TF/(sqrt(2)))^2/R_l;
P_ave_diff = diff(P_ave_output);
P_ave_diff = simplify(P_ave_diff, 'Steps',500)                      % <— Add This
P_ave_max = P_ave_diff == 0
R_l_optimal = solve(P_ave_max, R_l)
figure
fplot(P_ave_output,[10 200])
xticks([10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200])
xline(double(R_l_optimal))
ylabel('Average Power (Watts)')
xlabel('Load Resistance R_L (Ohms)')
figure
fplot(P_ave_diff,[1 2000])
ylabel('d(Average Power (Watts))/d(R_L)')
xlabel('Load Resistance R_L (Ohms)')
.
Related Question