Assigning values for a variable and solving symbolic equations numerically

Hello, how are you?

I am trying to assign values for a variable 't' to solve a symbolic equation. I used the command solve to solve the equations algebrically to the desired variable 'd0' and I would like to assign variables for t and obtain 'd0' numerically. How can I do this ?

Please find my code below:

syms d0 t positive
h = 10; % Height of the tank [m]
H = 30; % Height of the base [m]
D = 10; % Diameter of the tank [m]
w = 700; % Wind pressure of the tank [N/m²]
e = 10/100; % Eccentricity [m]
delta_a = 20/100; % Allowable deflection [m]
gamma_w = 10*10^3; % Unit weight of water [N/m³]
gamma_s = 80*10^3; % Unit weight of steel [N/m³]
E = 210*10^9; % Modulus of elasticity [Pa]
t_t = 1.5/100; % Average thickness of the tank wall
sigma_b = 165 * 10^6; % Allowable bending stress [Pa]
R = d0/2 - t;
I = pi/64*(d0^4 - (d0 - 2*t)^4);
A = pi*t*(d0 - t);
r = sqrt(I/A); % radius of gyration [m]
sigma_a = (12*pi^2*E)/(92*(H/r)^2); % Allowable axial stress [Pa]
V = 1.2*pi*D^2 * h; % Volume of the tank [m³]
A_s = 1.25*pi*D^2; % Surface area of tank [m²]
A_p = (2*D*h)/3; % Projected area for wind loading [m²]
P = V*gamma_w + A_s*t_t*gamma_s; % Load on the column due to weight of water and steel tank
W = w*A_p; % Lateral load at the tank C.G [m]
delta1 = (W*H^2)/(12*E*I) * (4*H + 3*h);
delta2 = H/(2*E*I) * (0.5*W*h + P*e)*(H + h);
delta = delta1 + delta2; % Deflection at C.G of the tank
M = W*(H + 0.5*h) + (delta + e)*P; % Moment at base [N.m]
f_b = M/(2*I) *d0; % Bending stress
f_a = V*gamma_w + A_s*t_t*gamma_s; % Axial stress
t = [1/100:0.01:20/100]; % t range #ok<NBRAK>
g1 = sigma_a - f_a == 0;
g1_d0 = solve(g1,d0)
g2 = sigma_b - f_b == 0;
g2_d0 = solve(g2,d0)
g3 = delta_a - delta ==0;
g3_d0 = solve(g3,d0)
g4 = R - 20 == 0;
g4_d0 = solve(g4,d0)
g5 = -R + 0.35 == 0;
g5_d0 = solve(g5,d0)

Thanks in advance !

Answer 1

[Paul]

You can use subs() and double(), at least for the "easy" versions of d0:

syms d0 t positive
h = 10; % Height of the tank [m]
H = 30; % Height of the base [m]
D = 10; % Diameter of the tank [m]
w = 700; % Wind pressure of the tank [N/m²]
e = 10/100; % Eccentricity [m]
delta_a = 20/100; % Allowable deflection [m]
gamma_w = 10*10^3; % Unit weight of water [N/m³]
gamma_s = 80*10^3; % Unit weight of steel [N/m³]
E = 210*10^9; % Modulus of elasticity [Pa]
t_t = 1.5/100; % Average thickness of the tank wall
sigma_b = 165 * 10^6; % Allowable bending stress [Pa]
R = d0/2 - t;
I = pi/64*(d0^4 - (d0 - 2*t)^4);
A = pi*t*(d0 - t);
r = sqrt(I/A); % radius of gyration [m]
sigma_a = (12*pi^2*E)/(92*(H/r)^2); % Allowable axial stress [Pa]
V = 1.2*pi*D^2 * h; % Volume of the tank [m³]
A_s = 1.25*pi*D^2; % Surface area of tank [m²]
A_p = (2*D*h)/3; % Projected area for wind loading [m²]
P = V*gamma_w + A_s*t_t*gamma_s; % Load on the column due to weight of water and steel tank
W = w*A_p; % Lateral load at the tank C.G [m]
delta1 = (W*H^2)/(12*E*I) * (4*H + 3*h);
delta2 = H/(2*E*I) * (0.5*W*h + P*e)*(H + h);
delta = delta1 + delta2; % Deflection at C.G of the tank
M = W*(H + 0.5*h) + (delta + e)*P; % Moment at base [N.m]
f_b = M/(2*I) *d0; % Bending stress
f_a = V*gamma_w + A_s*t_t*gamma_s; % Axial stress
tvals = [1/100:0.01:20/100]; % t range #ok<NBRAK>  % so as not to conflict with the sym variable t
g1 = sigma_a - f_a == 0;
g1_d0 = solve(g1,d0,'ReturnConditions',true)

g1_d0 = struct with fields:
            d0: [2×1 sym]
    parameters: [1×0 sym]
    conditions: [2×1 sym]

g2 = sigma_b - f_b == 0;
g2_d0 = solve(g2,d0,'ReturnConditions',true)

g2_d0 = struct with fields:
            d0: z1
    parameters: z1
    conditions: 260082878480646144000000000000000*pi^2*t^8 + 1820580149364523008000000000000000*pi^2*t^6*z1^2 + 1105352233542746112000000000000000*pi^2*t^4*z1^4 + 65020719620161536000000000000000*pi^2*t^2*z1^6 + 17182393990160644319477760000000*pi*t^4…

g3 = delta_a - delta ==0;
g3_d0 = solve(g3,d0,'ReturnConditions',true)

g3_d0 = struct with fields:
            d0: x
    parameters: x
    conditions: - 112742891520000000*pi*t^4 + 225485783040000000*pi*t^3*x - 169114337280000000*pi*t^2*x^2 + 56371445760000000*pi*t*x^3 - 31731444346091383 == 0 & 0 < x

g4 = R - 20 == 0;
g4_d0 = solve(g4,d0,'ReturnConditions',true)

g4_d0 = struct with fields:
            d0: 2*t + 40
    parameters: [1×0 sym]
    conditions: symtrue

g5 = -R + 0.35 == 0;
g5_d0 = solve(g5,d0,'ReturnConditions',true)

g5_d0 = struct with fields:
            d0: 2*t + 7/10
    parameters: [1×0 sym]
    conditions: symtrue

double(subs(g5_d0.d0,t,sym(tvals)))

ans = 1×20

0.7200 0.7400 0.7600 0.7800 0.8000 0.8200 0.8400 0.8600 0.8800 0.9000 0.9200 0.9400 0.9600 0.9800 1.0000 1.0200 1.0400 1.0600 1.0800 1.1000

You'll have to decide what you want to do with g1_d0, which has two solutions, and g2_d0 and g3_d0 which have to satisfy conditions.

[Rodrigo Pena]

Hello Paul, thanks for your answer !

I did not understant the 'returnConditions' and the struct with fields. Do I have to simply compy and paste the lines to my code as it is ?

Regarding g1_d0 I want to discard the first solution and keep the second, I took a look on documentation about returnconditions but I could not understand ):

[Paul]

You can access the second solution of g1_d0 by

g1_d0.d0(2)

Some equations only have solutions when certain conditions are true. Consdier the equation

(x+y)/(x+z) = 0.

This equation has a trivial solution, but only as long as y is not the same as z, which is what solve() show

syms x y z
sol = solve((x+y)/(x+z)==0,'ReturnConditions',true)

sol = struct with fields:
             x: -y
    parameters: [1×0 sym]
    conditions: y ~= z

Other equations have more than one (or infinite) solutions that can still be expressed in terms of free parameter. Consider the equation

sin(x) = 0

Any value of x that is an integer multiple of pi is solution, which is what solve() shows

sol = solve(sin(x)==0,'ReturnConditions',true)

sol = struct with fields:
             x: pi*k
    parameters: k
    conditions: in(k, 'integer')

Note here that sol.conditions shows the restriction that the parameter, k, must be an integer.

I didn't look too closely at the equations to see what's going on with g2_d0 and g3_d0.

[Rodrigo Pena]

aah ok. I managed to implement what you wrote above. I will search more on how to get rid of these z1 and x variables that appears in the calculations. thanks !