empty sym 0-by-1

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)')

Answer 1

[Star Strider]

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)')

.

[Shane Palmer]

Wow, Star Strider, what does that line of code do? Thanks for the fix. I noticed that it output my R_L_optimal in a complex number so I had to add an abs() command to reach the same R_L value.

Thanks a ton.

Shane

[Star Strider]

As always, my pleasure!

The solve function apparently has problems parsing some more complicated expressions, and so will just give up. (I haven’t experimented with Maple or any of the other symbolic engines.) Simplifying the expression (here telling it to keep simplifying until it can’t simplify further or reaches the iteration limit) solves that by giving solve a more tractable expression. Then, it works.

Interesting about ‘R_L_optimal’. When I ran it, I got:

R_l_optimal =
72.275099696632748765256742994398

a real number. (I suspect the imaginary part is vanishingly small.) This with R2020a, Update 2.