# Assigning values for a variable and solving symbolic equations numerically

6 views (last 30 days)
Rodrigo Pena on 13 Oct 2021
Commented: Rodrigo Pena on 14 Oct 2021
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 ?
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²]
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)

Paul on 13 Oct 2021
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²]
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 on 14 Oct 2021
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 !

R2020a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!