Do I have a convergence issue?

61 次查看（过去 30 天）
Timothy
Timothy 2024-11-11，9:02
回答： Torsten 2024-11-11，11:10
My code does not seem to find the correct B.
Can someone tell me where the code is wrong?
g = 10; %m/s^2
Q = 70*1000; %N/m
tau_L1 = 24*1000; %Pa
t_w = 300/1000; %m
t_f = 300/1000; %m
h_w = 2200; %mm
H_w = h_w/1000; %m
H_tot = H_w+t_f; %m
D_f = t_f; %m
rho_dim = 2.1*1000; %kg/m^3
rho_c = 2.3*1000; %kg/m^3
rho_L1 = 1.8*1000; % kg/m^3
p_max = rho_dim*g*H_tot; % Pa
phi_deg = 30; %degrees
phi_prim = deg2rad(phi_deg); %radians
K_0 = (1-sin(phi_prim))*2^(sin(phi_prim)); %no unit
P = p_max*(H_tot/2)*K_0; %N/m
G_c1 = t_w*H_w*rho_c*g; % N/m
B = 0.5; %m
tolerance = 1e-6; % no unit
max_iter = 3000; %no unit
iter = 0; % no unit
while iter < max_iter
iter = iter+1; %no unit
phi = (0.2*phi_prim+0*(B-0.2))/B; %radians
G_c2 = t_w*B*rho_c*g; %N/m
G_L = ((B-t_w)/2)*H_w*rho_dim*g; %N/m
alpha = atan(P/(Q+G_c1+G_c2+G_L)); %radians
x = (P*((t_w/3)+0.3)+G_L*((B-t_w)/2)+Q*(B/2)+G_c1*(B/2)+G_c2*(B/2))/(Q+G_c1+G_c2+G_L); %m
e = x-(B/2); %m
c_prim_grus = 0; %Pa
c_prim = (0.2*c_prim_grus+(B-0.2)*tau_L1)/B; %Pa
B_prim = B-2*e; %m
N_q = exp(pi*tan(phi))*tan((deg2rad(45))+((phi)/2))^2; %no unit
N_c = (N_q-1)*cot(phi); % no unit
N_gamma = 2*(N_q+1)*tan(phi); % no unit
lambda_ci = 1-(rad2deg(alpha)/90)^2; % no unit
lambda_qi = 1-(rad2deg(alpha)/90)^2; % no unit
lambda_gammai = (1-(alpha/phi))^2; % no unit
gamma = (g*rho_dim*0.2+g*(rho_L1-(1*10^3))*(B-0.2))/B; %Pa
R = sqrt(P^2+(Q+G_c1+G_c2+G_L)^2); %N/m
q = rho_dim*g*D_f; %N/m
if D_f/B_prim < 1
lambda_cd1 = 1 + 0.4*(D_f/B_prim); %no unit
lambda_gammad1 = 1; %no unit
lambda_qd1 = 1+2*tan(phi)*((1-sin(phi))^2)*(D_f/B_prim); %no unit
else
lambda_cd1 = 1 + 0.4*(D_f/B_prim); %no unit
lambda_gammad1 = 1; %no unit
lambda_qd1 = 1+2*tan(phi)*((1-sin(phi))^2)*atan(deg2rad(D_f/B_prim)); %no unit
end
lhs = (B-2*(e))/cos(alpha)*(c_prim*lambda_cd1*lambda_ci*N_c+q*lambda_qd1*lambda_qi*N_q+0.5*lambda_gammad1*lambda_gammai*gamma*(B-2*e)*N_gamma); %N/m
rhs = R; %N/m
residual = lhs-rhs; %N/m
if abs(residual) < tolerance
disp('Convergence reached')
break
end
B = B + 0.001; %m
end
disp(B)

请先登录，再回答此问题。

采纳的回答

Epsilon
Epsilon 2024-11-11，11:02
Hi Timothy,
The code is not converging and completing the max iteration steps(3000). This is happening as the "residual (lhs - rhs)" is of much higher order and failing the check “abs(residual) < tolerance” if tolerance is 1e-6.
To solve it there can be two approaches:
1. Keep the tolerance at higher order, around 100(<0.1% of LHS)
2. Increase the "max_iter" to a higher value and decrease the increment of "B". This will give a more accurate value of B but will take more time to calculate.
Result when values at max_iter=30000000, tolerance=1e-2 and increment of B by B = B + 0.0000001
Hope it helps.

更多回答（1 个）

Torsten
Torsten 2024-11-11，11:10
Ran in:
B0 = 0.5;
B = fsolve(@fun,B0)
Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.
B = 0.9247
function residual = fun(B)
g = 10; %m/s^2
Q = 70*1000; %N/m
tau_L1 = 24*1000; %Pa
t_w = 300/1000; %m
t_f = 300/1000; %m
h_w = 2200; %mm
H_w = h_w/1000; %m
H_tot = H_w+t_f; %m
D_f = t_f; %m
rho_dim = 2.1*1000; %kg/m^3
rho_c = 2.3*1000; %kg/m^3
rho_L1 = 1.8*1000; % kg/m^3
p_max = rho_dim*g*H_tot; % Pa
phi_deg = 30; %degrees
phi_prim = deg2rad(phi_deg); %radians
K_0 = (1-sin(phi_prim))*2^(sin(phi_prim)); %no unit
P = p_max*(H_tot/2)*K_0; %N/m
G_c1 = t_w*H_w*rho_c*g; % N/m
%B = 0.5; %m
%tolerance = 1e-6; % no unit
%max_iter = 3000; %no unit
%iter = 0; % no unit
%while iter < max_iter
% iter = iter+1; %no unit
phi = (0.2*phi_prim+0*(B-0.2))/B; %radians
G_c2 = t_w*B*rho_c*g; %N/m
G_L = ((B-t_w)/2)*H_w*rho_dim*g; %N/m
alpha = atan(P/(Q+G_c1+G_c2+G_L)); %radians
x = (P*((t_w/3)+0.3)+G_L*((B-t_w)/2)+Q*(B/2)+G_c1*(B/2)+G_c2*(B/2))/(Q+G_c1+G_c2+G_L); %m
e = x-(B/2); %m
c_prim_grus = 0; %Pa
c_prim = (0.2*c_prim_grus+(B-0.2)*tau_L1)/B; %Pa
B_prim = B-2*e; %m
N_q = exp(pi*tan(phi))*tan((deg2rad(45))+((phi)/2))^2; %no unit
N_c = (N_q-1)*cot(phi); % no unit
N_gamma = 2*(N_q+1)*tan(phi); % no unit
lambda_ci = 1-(rad2deg(alpha)/90)^2; % no unit
lambda_qi = 1-(rad2deg(alpha)/90)^2; % no unit
lambda_gammai = (1-(alpha/phi))^2; % no unit
gamma = (g*rho_dim*0.2+g*(rho_L1-(1*10^3))*(B-0.2))/B; %Pa
R = sqrt(P^2+(Q+G_c1+G_c2+G_L)^2); %N/m
q = rho_dim*g*D_f; %N/m
if D_f/B_prim < 1
lambda_cd1 = 1 + 0.4*(D_f/B_prim); %no unit
lambda_gammad1 = 1; %no unit
lambda_qd1 = 1+2*tan(phi)*((1-sin(phi))^2)*(D_f/B_prim); %no unit
else
lambda_cd1 = 1 + 0.4*(D_f/B_prim); %no unit
lambda_gammad1 = 1; %no unit
lambda_qd1 = 1+2*tan(phi)*((1-sin(phi))^2)*atan(deg2rad(D_f/B_prim)); %no unit
end
lhs = (B-2*(e))/cos(alpha)*(c_prim*lambda_cd1*lambda_ci*N_c+q*lambda_qd1*lambda_qi*N_q+0.5*lambda_gammad1*lambda_gammai*gamma*(B-2*e)*N_gamma); %N/m
rhs = R; %N/m
residual = lhs-rhs; %N/m
% if abs(residual) < tolerance
% disp('Convergence reached')
% break
% end
% B = B + 0.001; %m
end

类别

在 Help Center 和 File Exchange 中查找有关 Legend 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by