%I encountered the following problem in the calculation: 1. The calculated H is negative, %and I am unsure if the calculation is correct. Some formulas cannot be simplified and still %exist in the form of 5000/51166. 3. Poor overall code fluency clea

clear all
close
%% 参数定义parameter definition
P = 42;
c = 800;
E = 15000000;
K = 1.8;
P_ya = 18000;
F = 2;
y = 26.8;
R = 20; % radius
syms H B
%H=100;
% 计算破裂角 a Calculate the rupture angle
if K <= 0.5
    a = 90;
elseif K <= 1
    a = -90 * K + 135;
elseif K <= 3
    a = -22.5 * K + 67.5;
else
    a = 0;
end

% 显示计算得到的 a 的值
disp(['当 K = ', num2str(K), ' 时，破裂角 a = ', num2str(a), '°']);

%% 求解初始破裂角相关量 Solving the initial rupture angle related quantities
L = H + R * (1 - sind(a));
G_1 = (y * L^2) / (2 * tand(B)); % 三角形块体的自重
p = atan2(tand(P), F); % 折减后的内摩擦角
p = rad2deg(p); % 将结果转换为度,便于后面统一计算
C = c * L; % 竖直面上粘聚力合力
C_s = c * L / (F * sind(B)); % 破裂面上粘聚力合力
G_0 = 2 * y * H * cosd(a);
z = 0.9 * P; % 按照围岩等级取值，三级围岩取0.9

%% 定义目标函数 E(B) Define the objective function E (B)
%E_func = @(B) (y ./ (2 .* tand(B))) .* sind(B + p) ./ cosd(B + p - z);
E_func=@(B) (cosd(B+p).*sind(B)).*cosd(B).*cosd(B+p-z)+sind(B+p).*sind(B).*(sind(B+p-z).*cosd(B)+cosd(B+p-z).*sind(B));
%% 数值求导函数 Numerical derivative function
% 使用中心差分法计算导数
dE_func = @(B) (E_func(B + 1e-6) - E_func(B - 1e-6)) / (2e-6);
%% 数值寻找导数为零的 B 值
% 只寻找一个接近的 B 值
B_range = [0, 90]; % B 的取值范围
B_init = 45;           % 初始猜测值，设置为 45 度
% 使用 fzero 寻找导数为零的 B 值
try
    B_zero = fzero(dE_func, B_init);
    % 检查找到的 B 值是否满足条件
    if abs(dE_func(B_zero)) < 1e-6
        disp(['找到满足条件的 B 值为：', num2str(B_zero)]);
    else
        disp('没有找到导数接近零的 B 值');
    end
catch
    disp('fzero 计算失败，未找到满足条件的 B 值');
end
B=B_zero 
%% 计算埋深 Calculate burial depth
f1 = ((G_1 - C) .* sind(p + B_zero) + C_s .* cosd(p)) ./ cosd(B_zero + p - z);
f2 = (P_ya - G_0 - 2 .* C) / (2 * sind(z));
% 定义控制方程，解出 H
eqn = f1 - f2 == 0;
% 使用 solve 反解出 H
sol_H_sym = solve(eqn, H);
% 将符号解转换为具体的数值
sol_H_num = double(subs(sol_H_sym));
% 显示结果
disp(['解出的 H 的值为：', num2str(sol_H_num)]);

[Image Analyst]

I never use symbolic. Are you absolutely sure you must have H and B be symbolic? Can you do it numerically instead?

[Nicholas]

I didn't understand what you meant, could you be more specific?This will be very helpful to me

Answer 1

[Shashi Kiran]

Hi Nicholas,

I understand that you are trying to solve for Burial Depth(H) using the Symbolic Math Toolbox in MATLAB but not getting the expected results.

Here is an alternative approach using numerical methods to solve for H highlighting the modifications made to your code.

• The initializations, the definition of the objective function, the numerical derivative and calculation of B remain unchanged.

clear;close all;clc
% Parameter definition
P = 42;
c = 800;
E = 15000000;
K = 1.8;
P_ya = 18000;
F = 2;
y = 26.8;
R = 20; % radius

% Calculate the rupture angle
if K <= 0.5
    a = 90;
elseif K <= 1
    a = -90 * K + 135;
elseif K <= 3
    a = -22.5 * K + 67.5;
else
    a = 0;
end

% Display the calculated value of a
disp(['When K = ', num2str(K), ', the rupture angle a = ', num2str(a), '°']);

When K = 1.8, the rupture angle a = 27°

p = atan2(tand(P), F); % Reduced internal friction angle
p = rad2deg(p); % Convert the result to degrees for uniform calculation later
z = 0.9 * P; % Value based on rock mass grade, 0.9 for grade 3 rock mass
E_func = @(B) (cosd(B+p).*sind(B)).*cosd(B).*cosd(B+p-z)+sind(B+p).*sind(B).*(sind(B+p-z).*cosd(B)+cosd(B+p-z).*sind(B));
dE_func = @(B) (E_func(B + 1e-6) - E_func(B - 1e-6)) / (2e-6);

B_init = 45; % Initial guess, set to 45 degrees

% Use fzero to find the B value where the derivative is zero
try
    B_zero = fzero(dE_func, B_init);
    disp(['Found B value that meets the condition: ', num2str(B_zero)]);
catch
    error('fzero computation failed, did not find a B value that meets the condition');
end

Found B value that meets the condition: 56.1005

• Defined a function f that represents the equation to solve for H. Used fsolve with an initial guess H_init.
• Removed the symbolic solution approach (solve and subs) and directly calculated H numerically.

% Define the function to find zero
f = @(H) ((y * (H + R * (1 - sind(a)))^2) / (2 * tand(B_zero)) - c * (H + R * (1 - sind(a)))) * sind(p + B_zero) + ...
         (c * (H + R * (1 - sind(a))) / (F * sind(B_zero))) * cosd(p) / cosd(B_zero + p - z) - ...
         (P_ya - 2 * y * H * cosd(a) - 2 * c) / (2 * sind(z));

% Initial guess for H
H_init = 0;

% Solve for H using fsolve
options = optimset('Display', 'off'); % Suppress output
H_solution = fsolve(f, H_init, options);

% Display the result
disp(['The calculated value of H is: ', num2str(H_solution)]);

The calculated value of H is: 38.0848

Refer to the following documentations for more details about the functions:

1. fsolve: https://www.mathworks.com/help/optim/ug/fsolve.html
2. function_handle: https://www.mathworks.com/help/matlab/ref/function_handle.html

Hope this helps.

[Nicholas]

G_1 = (y * L^2) / (2 * tand(B));

E_func = @(B) (y.*L^2 ./ (2 .* tand(B))) .* sind(B + p) ./ cosd(B + p - z);

I got confused. This is Ea, and above it is the derivative dE.

I need to first take the derivative of E0 to obtain B0, then use formula 13 to obtain H1, and then use formula 16 to obtain the next B1. Compare the difference between B0 and B1 and iterate repeatedly

thanks