How do I make solve() return solutions that aren't empty?

syms D d 
assume(d > 0 & D > d)
for n = 1:length(omega)
    I = (pi/64)*(D^4 - d^4);
    k = (3*E*I) / (l^3);
    c = zeta*2*sqrt(m*k);
    z = F ./ (k-m.*omega.^2) + (c.^2).*(omega.^2).^0.5;
    sigma_1 = (3*D*z) / (2*l^2);
    sigma_2 = (m*g) / ((pi/4) * (D^2 - d^2));
    M_result = (pi/4) * (D^2 - d^2)*l*rho_steel; 
    omega_n = sqrt(k/m);
    
    eq1 = sigma_1 + sigma_2 == E;
    eq2 = omega_n == omega_n_min;
    eq3 = M_result == M;
    
    solution = solve([eq1, eq2, eq3], D, d);
    solution.D
    solution.d
end

I have been searching through other threads with people asking similar questions to mine, but I can't figure out why I am only getting empty solutions. I am not getting any errors, but my solutions for 'D' and 'd' are still empty. What should I change to fix this issue?

[Stephan]

Do you have values for

E, l, zeta, m, F, omega, g, rho_steel

?

[madhan ravi]

solving using solve inside loop is inefficient

[Stephen Doroba]

Yes I do have values for all of those. I tried doing it outside of a for loop as well, but the answer is empty regardless.

[Stephan]

See Walters answer - you forgot the indices of omega when running the loop. This is problably the solution to your problem:

z = F ./ (k-m.*omega(n).^2) + (c.^2).*(omega(n).^2).^0.5;

[Stephen Doroba]

I made those changes and the result is the same.

syms D d
assume(d > 0 & D > d)
for n = 1:length(omega)
    E = E_given * 6894.76; % convert to Pa
    I = (pi/64)*(D^4 - d^4);
    k = (3*E*I) / (l^3);
    c = zeta*2*sqrt(m*k);
    z = F ./ ((k-m.*omega(n).^2) + (c.^2)*(omega(n).^2))^0.5;
    sigma_1 = (3*D*z) / (2*l^2);
    sigma_2 = (m*g) / ((pi/4) * (D^2 - d^2));
    M_result = (pi/4) * (D^2 - d^2)*l*rho_steel; 
    omega_n = sqrt(k/m);
    
    eq1 = sigma_1 + sigma_2 == E;
    eq2 = omega_n == omega_n_min;
    eq3 = M_result == M;
    
    solution = solve([eq1, eq2, eq3], [D, d]);
end

solution.D
solution.d
   

ans =
Empty sym: 0-by-1
ans =
Empty sym: 0-by-1
>> 

[Stephan]

could you provide the missing values?

[Stephen Doroba]

Apologies, here is the full code:

clc, clear all, format compact;

W_given = 100000; % lb
l_given = 50; % ft
E_given = 30000; % psi
g = 9.81; % m/s

A = 0.5*g;
omega = 2*pi*[5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20];

m = W_given * 0.453592; % convert to kg
l = l_given * 0.3048; % convert to m
E = E_given * 6894.76; % convert to Pa
rho_steel = 8050; % kg/m^3
F = A*m; % N 
zeta = 0.15;
M = 200000; % Mass of support (kg)
omega_n_min = 20*2*pi;

syms D d
assume(d > 0 & D > d)
for n = 1:length(omega)
    E = E_given * 6894.76; % convert to Pa
    I = (pi/64)*(D^4 - d^4);
    k = (3*E*I) / (l^3);
    c = zeta*2*sqrt(m*k);
    z = F ./ ((k-m.*omega(n).^2) + (c.^2)*(omega(n).^2))^0.5;
    sigma_1 = (3*D*z) / (2*l^2);
    sigma_2 = (m*g) / ((pi/4) * (D^2 - d^2));
    M_result = (pi/4) * (D^2 - d^2)*l*rho_steel; 
    omega_n = sqrt(k/m);
    
    eq1 = sigma_1 + sigma_2 == E;
    eq2 = omega_n == omega_n_min;
    eq3 = M_result == M;
    
    solution = solve([eq1, eq2, eq3], [D, d]);
end

diameters = single(solution.D)

Answer 1

[Walter Roberson]

when you construct z you need omega(n)

[Stephen Doroba]

That results in the same issue.

syms D d

for n = 1:length(omega)
    E = E_given * 6894.76; % convert to Pa
    I = (pi/64)*(D^4 - d^4);
    k = (3*E*I) / (l^3);
    c = zeta*2*sqrt(m*k);
    z = F / ((k-m*omega(n)^2) + (c^2)*(omega(n)^2))^0.5;
    sigma_1 = (3*D*z) / (2*l^2);
    sigma_2 = (m*g) / ((pi/4) * (D^2 - d^2));
    M_result = (pi/4) * (D^2 - d^2)*l*rho_steel; 
    omega_n = sqrt(k/m);
    
    eq1 = sigma_1 + sigma_2 == E;
    eq2 = omega_n == omega_n_min;
    eq3 = M_result == M;
    
    solution = solve([eq1, eq2, eq3], [D, d]);
end

solution.D
solution.d

ans =
Empty sym: 0-by-1
ans =
Empty sym: 0-by-1

[madhan ravi]

so provide the values for us to experiment

[Stephen Doroba]

Apologies, here is everything:

clc, clear all, format compact;

W_given = 100000; % lb
l_given = 50; % ft
E_given = 30000; % psi
g = 9.81; % m/s

A = 0.5*g;
omega = 2*pi*[5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20];

m = W_given * 0.453592; % convert to kg
l = l_given * 0.3048; % convert to m
E = E_given * 6894.76; % convert to Pa
rho_steel = 8050; % kg/m^3
F = A*m; % N 
zeta = 0.15;
M = 200000; % Mass of support (kg)
omega_n_min = 20*2*pi;

syms D d
assume(d > 0 & D > d)
for n = 1:length(omega)
    E = E_given * 6894.76; % convert to Pa
    I = (pi/64)*(D^4 - d^4);
    k = (3*E*I) / (l^3);
    c = zeta*2*sqrt(m*k);
    z = F ./ ((k-m.*omega(n).^2) + (c.^2)*(omega(n).^2))^0.5;
    sigma_1 = (3*D*z) / (2*l^2);
    sigma_2 = (m*g) / ((pi/4) * (D^2 - d^2));
    M_result = (pi/4) * (D^2 - d^2)*l*rho_steel; 
    omega_n = sqrt(k/m);
    
    eq1 = sigma_1 + sigma_2 == E;
    eq2 = omega_n == omega_n_min;
    eq3 = M_result == M;
    
    solution = solve([eq1, eq2, eq3], [D, d]);
end

diameters = single(solution.D)
   

Answer 2

[Stephan]

Hi,

i think your problem does not have a feasible solution. If you rewrite it a little bit and use fsolve the resulting message is "No solution found". This is the code i used to find a numerical solution to your problem:

result = solve_D_d

function res = solve_D_d
W_given = 100000; % lb
l_given = 50; % ft
E_given = 30000; % psi
g = 9.81; % m/s
A = 0.5*g;
omega = 2*pi*[5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20];
m = W_given * 0.453592; % convert to kg
l = l_given * 0.3048; % convert to m
E = E_given * 6894.76; % convert to Pa
rho_steel = 8050; % kg/m^3
F = A*m; % N
zeta = 0.15;
M = 200000; % Mass of support (kg)
omega_n_min = 20*2*pi;
res=zeros(2,numel(omega));
options = optimoptions(@fsolve,'Algorithm','Levenberg-Marquardt','StepTolerance',1e-12)
for i = 1:numel(omega)
    res(:,i) = fsolve(@calc_D_d,[2 1],options);
end
calc_D_d(res(:,1))
    function Fun = calc_D_d(x)
    D = x(1);
    d = x(2);
    
    k = (3*E*(pi/64)*(D^4 - d^4)) / (l^3);
    c = zeta*2*sqrt(m*k);
    z = F ./ ((k-m.*omega(i).^2) + (c.^2)*(omega(i).^2))^0.5;
    sigma_1 = (3*D*z) / (2*l^2);
    sigma_2 = (m*g) / ((pi/4) * (D^2 - d^2));
    M_result = (pi/4) * (D^2 - d^2)*l*rho_steel;
    omega_n = sqrt(k/m);
    
    Fun(1) = sigma_1 + sigma_2 - E;
    Fun(2) = omega_n - omega_n_min;
    Fun(3) = M_result - M;
end  
end

Unfortunately this is also not successful and does not find a viable solution.

What could you do now?

1. You could try the suggestions given in the documentation if fsolve fails. But i guess that there is indeed no solution to the problem, due to the way you formulated it.
2. Another approach would be to use fgoalattain combined with restrictions for upper and lower goals. IMO the problem you have is that you want to meet 3 goals exactly, which appears to be not possible. So you could reformulate the problem by trying to meet a target range of your values for the 3 functions you defined. For example you could define that sigma1 + sigma2 >= 0.9*E and sigma1 + sigma2 <= 1.1*E. I could imagine that you will have more luck with this kind of approach. An example how to define this using fgoalattain can be found here.
3. Maybe other optimizers are also suitable to find a feasible solution inside a target range of your functions.

Best regards

Stephan