Do I have a convergence issue?

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)

Answer 1

[Epsilon]

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.

Answer 2

[Torsten]

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