用二分法求解A*sin(w*x)+k*exp(L*x)代码

时间: 2024-06-22 14:02:29 浏览: 126
使用二分法(Binary Search)求解非线性方程 `A*sin(w*x) + k*exp(L*x) = 0` 的数学问题时,你需要编写一个迭代算法,该算法在给定一定精度要求的情况下,不断缩小搜索范围,直到找到方程的根。这是一种数值方法,通常用于没有封闭形式解的问题。 以下是 Python 中的一个基本示例,使用 `scipy.optimize.bisect` 函数实现二分法: ```python import numpy as np from scipy.optimize import bisect def equation(x, A, w, k, L): return A * np.sin(w * x) + k * np.exp(L * x) def bisection_method(A, w, k, L, initial_guess, tolerance=1e-6, max_iterations=100): def interval_search(start, end): guess = start + (end - start) / 2 if abs(equation(guess, A, w, k, L)) < tolerance: return guess elif equation(start, A, w, k, L) * equation(guess, A, w, k, L) > 0: # 如果在开始和猜测值之间,向左搜索 return interval_search(start, guess) else: # 否则,在猜测值和结束之间搜索 return interval_search(guess, end) root = interval_search(initial_guess, 1e6) # 初始区间设置为足够大的数防止溢出 if abs(equation(root, A, w, k, L)) > tolerance: print(f"未能找到满足条件的根,可能需要更精确的初始猜测或增加迭代次数") return root # 示例用法 A, w, k, L = 1, 2, 3, 4 # 假设这些是你的系数 initial_guess = 0.1 # 初始猜测的根位置 root = bisection_method(A, w, k, L, initial_guess) ```
阅读全文

相关推荐

import numpy as np import matplotlib.pyplot as plt # 材料参数假设(用户未提供) k = 0.5 # 导热系数 (W/m·K) rho = 1000 # 密度 (kg/m³) cp = 4000 # 比热容 (J/kg·K) alpha = k / (rho * cp) # 热扩散率 (m²/s) # 几何参数 Lx = 0.05 # x方向半长 (m) Ly = 0.01 # y方向半宽 (m) Lz = 0.02825 # z方向半高 (m) # 热物性参数 h = 500 # 对流换热系数 (W/m²·K) T0 = 20 # 初始温度 (°C) Tf = 100 # 介质温度 (°C) # 计算毕渥数 Bi_x = h * Lx / k Bi_y = h * Ly / k Bi_z = h * Lz / k # 时间数组 (0-60秒) time = np.linspace(0, 60, 61) def find_root(bi, n): """ 二分法求解特征方程 λ*tan(λ) = Bi """ a = (n - 1) * np.pi b = (n - 0.5) * np.pi - 1e-9 # 避开奇点 for _ in range(100): c = (a + b) / 2 fc = c * np.tan(c) - bi if fc < 0: a = c else: b = c if (b - a) < 1e-6: return (a + b) / 2 return (a + b) / 2 def calculate_theta(roots, L, t, position='center'): """ 计算一维无量纲温度 """ theta = 0.0 for lam in roots: # 系数项计算 C = (2 * np.sin(lam)) / (lam + np.sin(lam)*np.cos(lam)) # 指数项 Fo = alpha * t / L**2 exp_term = np.exp(-lam**2 * Fo) # 空间项 if position == 'center': space_term = 1.0 else: space_term = np.cos(lam) theta += C * exp_term * space_term # 收敛判断 if abs(C * exp_term * space_term) < 1e-4: break return theta # 求每个方向的前5个根 roots_x = [find_root(Bi_x, n) for n in range(1, 6)] roots_y = [find_root(Bi_y, n) for n in range(1, 6)] roots_z = [find_root(Bi_z, n) for n in range(1, 6)] # 计算温度曲线 Tc = [] # 中心温度 Ta = [] # 表面温度(假设三个方向表面) for t in time: # 三维中心温度 theta_x = calculate_theta(roots_x, Lx, t, 'center') theta_y = calculate_theta(roots_y, Ly, t, 'center') theta_z = calculate_theta(roots_z, Lz, t, 'center') Tc.append(Tf + (T0 - Tf) * theta_x * theta_y * theta_z) # 三维表面温度(三个方向表面) theta_x_s = calculate_theta(roots_x, Lx, t, 'surface') theta_y_s = calculate_theta(roots_y, Ly, t, 'surface') theta_z_s = calculate_theta(roots_z, Lz, t, 'surface') Ta.append(Tf + (T0 - Tf) * theta_x_s * theta_y_s * theta_z_s) # 绘图 plt.figure(figsize=(10, 6)) plt.plot(time, Tc, label='中心温度 Tc') plt.plot(time, Ta, label='表面温度 Ta', linestyle='--') plt.xlabel('时间 (秒)') plt.ylabel('温度 (℃)') plt.title('温度随时间变化曲线') plt.legend() plt.grid(True) plt.show() 帮我把以上代码转换成matlab语言

% 材料参数假设(用户未提供) k = 0.5; % 导热系数 (W/m·K) rho = 1000; % 密度 (kg/m³) cp = 4000; % 比热容 (J/kg·K) alpha = k / (rho * cp); % 热扩散率 (m²/s) % 几何参数 Lx = 0.05; % x方向半长 (m) Ly = 0.01; % y方向半宽 (m) Lz = 0.02825; % z方向半高 (m) % 热物性参数 h = 500; % 对流换热系数 (W/m²·K) T0 = 20; % 初始温度 (°C) Tf = 100; % 介质温度 (°C) % 计算毕渥数 Bi_x = h * Lx / k; Bi_y = h * Ly / k; Bi_z = h * Lz / k; % 时间数组 (0-60秒) time = linspace(0, 60, 61); % 求每个方向的前5个根 roots_x = arrayfun(@(n) find_root(Bi_x, n), 1:5); roots_y = arrayfun(@(n) find_root(Bi_y, n), 1:5); roots_z = arrayfun(@(n) find_root(Bi_z, n), 1:5); % 预分配数组 Tc = zeros(size(time)); % 中心温度 Ta = zeros(size(time)); % 表面温度(假设三个方向表面) % 计算温度曲线 for i = 1:length(time) t = time(i); % 三维中心温度 theta_x = calculate_theta(roots_x, Lx, t, 'center'); theta_y = calculate_theta(roots_y, Ly, t, 'center'); theta_z = calculate_theta(roots_z, Lz, t, 'center'); Tc(i) = Tf + (T0 - Tf) * theta_x * theta_y * theta_z; % 三维表面温度 theta_x_s = calculate_theta(roots_x, Lx, t, 'surface'); theta_y_s = calculate_theta(roots_y, Ly, t, 'surface'); theta_z_s = calculate_theta(roots_z, Lz, t, 'surface'); Ta(i) = Tf + (T0 - Tf) * theta_x_s * theta_y_s * theta_z_s; end % 绘图 figure('Position', [100 100 1000 600]) plot(time, Tc, 'LineWidth', 1.5) hold on plot(time, Ta, '--', 'LineWidth', 1.5) xlabel('时间 (秒)') ylabel('温度 (℃)') title('温度随时间变化曲线') legend('中心温度 Tc', '表面温度 Ta', 'Location', 'best') grid on hold off % ========== 自定义函数 ========== function lam = find_root(bi, n) % 二分法求解特征方程 λ*tan(λ) = Bi a = (n - 1) * pi; b = (n - 0.5) * pi - 1e-9; % 避开奇点 for iter = 1:100 c = (a + b)/2; fc = c * tan(c) - bi; if fc < 0 a = c; else b = c; end if (b - a) < 1e-6 lam = (a + b)/2; return end end lam = (a + b)/2; end function theta = calculate_theta(roots, L, t, position) % 计算一维无量纲温度 theta = 0.0; for j = 1:length(roots) lam = roots(j); % 系数项计算 C = (2 * sin(lam)) / (lam + sin(lam)*cos(lam)); % 指数项 Fo = alpha * t / L^2; exp_term = exp(-lam^2 * Fo); % 空间项 if strcmp(position, 'center') space_term = 1.0; else space_term = cos(lam); end theta = theta + C * exp_term * space_term; % 收敛判断 if abs(C * exp_term * space_term) < 1e-4 break end end end 以上是完整matlab代码,错误使用 alpha 输出参数太多。 出错 cft>calculate_theta (第 92 行) Fo = alpha * t / L^2; 出错 cft (第 39 行) theta_x = calculate_theta(roots_x, Lx, t, 'center');根据提示的修改意见,告诉我该怎么修改

clc;clear;close all; D=6; B=pi*D/4; gamma2=20; cc1_value=15; beita=0; C=2*D; fai_range=25:1:45; fai_value=deg2rad(fai_range); FOS_results = zeros(size(fai_value)); for ii=1:length(fai_value) fai=fai_value(ii); FOS0=1.0; toler=1e-3;%收敛误差 max_iter=1000;%最大迭代次数 step=0.1;%步长 for iter=1:max_iter cc1=cc1_value./FOS; fai1=atan(tan(fai)./FOS); delta=0.8*fai1; alpha2=pi/4+fai1/2; L=D/tan(alpha2); Ka=(1-sin(fai1))./(1+sin(fai1)); theta2=pi/4+fai1/2; f_11=(Ka-(Ka.*sin(theta2)^2)./3+sin(theta2)^2./3)./(sin(theta2)^2+Ka.*cos(theta2)^2); f_12=f_11.*cos(theta2)^2+sin(theta2)^2./3-1; f_31=(1-Ka).*cos(fai1)./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); f_32=cos(fai1)./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); f_21=(Ka.*(1+sin(fai1))+(1-sin(fai1)))./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); f_22=(cos(theta2)^2-sin(theta2)^2-sin(fai1))./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); Kanw0=1-sin(fai1) C01_D=C./D; Cpie=L/Kanw0; Ve=pi*B*Cpie*L/4; Se=B*L; C01=min(Ve./Se,C);%筒仓高度 P=2.*cc1.*sqrt(Ka).*(f_32-f_31.*f_12./f_11); E=f_31./f_11; J=2*(B+L)./(B*L); Q1=0+gamma2.*(C-C01); Q=(gamma2-J.*P)./(J.*E)+(Q1-(gamma2-J.*P)./(J.*E)).*exp(-J.*E.*C01);%筒仓传递下来的 theta=pi/2-1/2.*(asin(sin(delta)./sin(fai1)))+delta./2; Kanw1=3.*(cos(theta).^2+Ka.*sin(theta).^2)./(3-(1-Ka).*cos(theta).^2); Kanw2=3.*(cos(theta2).^2+Ka.*sin(theta2).^2)./(3-(1-Ka).*cos(theta2).^2); A1=Kanw1.*tan(delta);%tao_aw %tao_s A2=3*(1-Ka)*cos(fai1)/(2*(3-(1-Ka)*cos(theta)^2)); %seigema_s A3=3*(1-Ka)*cos(fai1)^2.*cot(alpha2)/(2*sin(fai1)*(3-(1-Ka)*cos(theta)^2)); %seigema_anw A4=Kanw1.*tan(beita); %tao_aw B1=cc1.*tan(delta).*(1+(1-Ka).*(((Kanw1*cos(theta).^2)/3)-sin(theta).^2))./tan(fai1);%tao_aw %tao_s B2=cc1*(1-Ka)*cos(fai1).^2./2*sin(fai1); B2=B2*(1+(1-Ka)*cos(theta)^2./3*(3-(1-Ka)*cos(theta)^2)); %seigema_s B3=cc1*(1-Ka)*cos(fai1).^2.*cot(alpha2)/(2*tan(fai1)*sin(fai1)); B3=B3*(1+(1-Ka)*cos(theta)^2./3*(3-(1-Ka)*cos(theta)^2))-cc1*cot(alpha2)./tan(fai1); %seigema_anw B4=cc1*(1-Ka)*tan(beita).*(((Kanw1*cos(theta).^2)/3)-sin(theta).^2)./tan(fai1); yita_1=A1+A2+A3+A4; yita_2=B1+B2+B3+B4; r=gamma2-2*cc1*sin(alpha2)./B;%对比时把后面的去掉 s=yita_1./(cot(alpha2)+tan(beita)); t=2*sin(alpha2)*Kanw2*tan(fai1)./B; kesi=yita_2./(cot(alpha2)+tan(beita)); tspan=[0 5.99]; f=@(z,q) [(r-kesi./(D-z))-((s-1)./(D-z)+t)*q] q0=Q; [z,q]=ode45(f, tspan, q0) seigema_h=Kanw1*q+cc1*(1-Ka).*(Kanw1*cos(theta).^2/3-sin(theta).^2)./tan(fai1);%水平土压力 if abs(seigema_h)<toler break; end grid off; box off; set(gca, 'LineWidth', 1.5, 'FontSize', 15, 'FontName', 'Times New Roman'); xlabel('$\varphi (^{\circ})$', 'Interpreter', 'latex', 'FontSize', 14, 'FontName', 'Times New Roman'); ylabel('$\sigma_{\mathrm{T}} / \gamma D$', 'FontSize', 14, 'FontName', 'Times New Roman', 'Interpreter', 'latex'); legend_labels = arrayfun(@(C) sprintf('$C/D=%.1f$', C/D), C_values, 'UniformOutput', false); legend(legend_labels, 'Location', 'northeast', 'Orientation', 'vertical', ... 'FontSize', 14, 'FontName', 'Times New Roman', 'Interpreter', 'latex', 'Box', 'off'); 这是我写的程序,可能有一些错误,帮我用简单迭代法求解安全系数FOS

clc;clear;close all; D=6; gamma2=20; B=pi*D/4; Q1=0;%地表荷载 fai_1=deg2rad(30); cc_1=2; cc1=@(FOS) cc_1./FOS; fai1=@(FOS) atan(tan(fai_1)./FOS); theta2=@(FOS) pi/4+fai1/2; alpha2=@(FOS) pi/4+fai1/2; L=@(FOS) D/tan(alpha2); Ka=@(FOS) (1-sin(fai1))./(1+sin(fai1)); f_11=@(FOS) (Ka-(Ka.*sin(theta2)^2)./3+sin(theta2)^2./3)./(sin(theta2)^2+Ka.*cos(theta2)^2); f_12=@(FOS) f_11.*cos(theta2)^2+sin(theta2)^2./3-1; f_31=@(FOS) (1-Ka).*cos(fai1)./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); f_32=@(FOS) cos(fai1)./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); f_21=@(FOS) (Ka.*(1+sin(fai1))+(1-sin(fai1)))./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); f_22=@(FOS) (cos(theta2)^2-sin(theta2)^2-sin(fai1))./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); Kanw0=@(FOS) f_21./f_11; C=10; Cpie=@(FOS) L/Kanw0; Ve=@(FOS) pi*B*Cpie*L/4; Se=@(FOS) B*L; C01=@(FOS) min(Ve./Se,C);%筒仓高度 beita=0; figure; hold on; P=@(FOS) 2.*cc1.*sqrt(Ka).*(f_32-f_31.*f_12./f_11); E=@(FOS) f_31./f_11; J=@(FOS) 2*(B+L)./(B*L); Q=@(FOS) (gamma2-J.*P)./(J.*E)+(Q1-(gamma2-J.*P)./(J.*E)).*exp(-J.*E.*C01); delta=@(FOS) 0.8*fai1; theta=@(FOS) pi/2-1/2.*(asin(sin(delta)./sin(fai1)))+delta./2; Kanw1=@(FOS) 3.*(cos(theta).^2+Ka.*sin(theta).^2)./(3-(1-Ka).*cos(theta).^2); Kanw2=@(FOS) 3.*(cos(theta2).^2+Ka.*sin(theta2).^2)./(3-(1-Ka).*cos(theta2).^2); A1=@(FOS) Kanw1.*tan(delta);%tao_aw %tao_s A2=@(FOS) 3*(1-Ka)*cos(fai1)/(2*(3-(1-Ka)*cos(theta)^2)); %seigema_s A3=@(FOS) 3*(1-Ka)*cos(fai1)^2.*cot(alpha2)/(2*sin(fai1)*(3-(1-Ka)*cos(theta)^2)); %seigema_anw A4=@(FOS) Kanw1.*tan(beita); %tao_aw B1=@(FOS) cc1.*tan(delta).*(1+(1-Ka).*(((Kanw1*cos(theta).^2)/3)-sin(theta).^2))./tan(fai1);%tao_aw %tao_s B2=@(FOS) cc1*(1-Ka)*cos(fai1).^2./2*sin(fai1); B2=@(FOS) B2*(1+(1-Ka)*cos(theta)^2./3*(3-(1-Ka)*cos(theta)^2)); %seigema_s B3=@(FOS) cc1*(1-Ka)*cos(fai1).^2.*cot(alpha2)/(2*tan(fai1)*sin(fai1)); B3=@(FOS) B3*(1+(1-Ka)*cos(theta)^2./3*(3-(1-Ka)*cos(theta)^2))-cc1*cot(alpha2)./tan(fai1); %seigema_anw B4=@(FOS) cc1*(1-Ka)*tan(beita).*(((Kanw1*cos(theta).^2)/3)-sin(theta).^2)./tan(fai1); yita_1=@(FOS) A1+A2+A3+A4; yita_2=@(FOS) B1+B2+B3+B4; r=@(FOS) gamma2-2*cc1*sin(alpha2)./B; s=@(FOS) yita_1./(cot(alpha2)+tan(beita)); t=@(FOS) 2*sin(alpha2)*Kanw2*tan(fai1)./B; kesi=@(FOS) yita_2./(cot(alpha2)+tan(beita)); tspan=[0 5.99]; f=@(z,q) [(r-kesi./(D-z))-((s-1)./(D-z)+t)*q] q0=Q; [z,q]=ode45(f, tspan, q0) seigama_anw=@(FOS) Kanw1*q+cc1*(1-Ka).*(Kanw1*cos(theta).^2/3-sin(theta).^2)./tan(fai1);%水平土压力 用cc1=@(FOS) cc_1./FOS; fai1=@(FOS) atan(tan(fai_1)./FOS);这两个求安全系数的公式,算一下,安全系数的曲线,具体思路是怎么算呢,是要循环安全系数FOS吗

% 椭圆轨道 clear all clc % 碎片参数 m0 = 10; % 碎片质量,单位 kg Cd = 2; % 阻力系数 A = 0.05; % 截面积,单位 m^2 B = m0 / (Cd * A); F_107 = 138; Ap = 4; % 轨道参数 Re = 6378e3; % 地球半径,单位 m hp = 247e3; % 近地点高度,单位 m ha = 247e3: 1e3: 2000e3; % 远地点高度,单位 m % rp0 = hp + Re; % 近地点半径,单位 m % ra = ha + Re; % 远地点半径,单位 m % a = (ra + rp0) / 2;% 轨道半长轴,单位 m % e = 1 - rp0 ./ a; % 偏心率,单位 rad mu = 3.986e14; % 引力参数m^3/s^2 f = 0; % 真近点角,单位 rad h_L = 133.8e3; % 再入高度,单位 m r_L = h_L + Re; % 再入半径,单位 m for i = length(ha) rp0 = hp + Re; % 近地点半径,单位 m ra0 = ha(i) + Re; % 远地点半径,单位 m a0 = (rp0 + ra0)/2; % 半长轴,单位 m e0 = 1 - rp0 / a0; % 偏心率 % 不同初始椭圆轨道的初始条件 y0 = [a0; e0]; tspan = [0 3000*24*3600];% 时间范围(初始猜测) % 事件函数:当r_p <= r_L时停止积分 options = odeset('Events', @(t,y) stopCondition1(t, y), 'RelTol', 1e-12,'MaxStep', 10000); [t, y] = ode45(@(t,y) dy_dt(t, y, B, mu, F_107, Ap), tspan, y0, options); % 记录碎片寿命 t_Life(i) = t(end) / (3600 * 24);% 单位 天 end % for i = 1: length(ha) % % rp0 = hp + Re; % 近地点半径,单位 m % ra0 = ha(i) + Re; % 远地点半径,单位 m % a0 = (rp0 + ra0)/2; % 半长轴,单位 m % e0 = 1 - rp0 / a0; % 偏心率 % % % 不同初始椭圆轨道的初始条件 % y0 = [0; a0; e0]; % 初始时间 % rp_span = [rp0 r_L]; % 近地点半径积分范围 % [rp, y] = ode45(@(rp,y) dt_dy(rp, y, B, mu, F_107, Ap), rp_span, y0); % % % 记录碎片寿命 % t_Life(i) = y(end,1); % % end plot(ha/1e3, t_Life); xlabel('远地点高度 (km)'); ylabel('轨道寿命 (天)');function dydt = dy_dt(t, y, B, mu, F_107, Ap) a = y(1); e = y(2); % rp = a * (1 - e); n = sqrt(mu / a^3); M0 = 0; t0 = 0; M = M0 + n * (t - t0); E = solve_kepler(M, e); f = 2 * atan2(sqrt((1+e)/(1-e)) * sin(E/2), cos(E/2)); % % 计算远地点半径 % r_a = a * (1 + e); % % % 检查并交换 r_p 和 r_a % if r_a < r_p % temp = r_p; % r_p = r_a; % r_a = temp; % end % % % 重新计算 a 和 e % a = (r_p + r_a) / 2; % e = 1 - r_p / a; v = sqrt(mu / a * (1 + 2 * e * cos(f) + e^2) / (1 - e^2)); % f = 0; % vp = sqrt(mu * (2 / r_p - 1 / a)); % 近地点轨道速度,单位 m/s h = a - 6378e3; % 大气模型 T = 900 + 2.5 * (F_107 - 70) + 1.5 * Ap; % 温度,单位 K m = 27 - 0.012 * (h / 1e3 - 200); % 摩尔质量 H = T *1000 / m; % 大气标高,单位 m p = 6e-10 * exp(-(h - 175e3) / H); % 大气密度,单位 kg/m^3 % p = 2.789e-10 * exp(-(h - 200e3) / 37.105e3); % 大气密度,单位 kg/m^3 % T0 = 2 * pi * sqrt(a^3 / mu);% 轨道周期 da_dt = a^2 * v^3 * p / (mu * B); de_dt = -v * (e + cos(f)) * p / B; % drp_dt = a * vp * p * (-1/mu * a * vp^2 * (1 - e) + (e + cos(f))) / B; % n = sqrt(mu / a^3); % da_dt = -(p / B) * n * a^2 * (1 + 2 * e + e^2)^1.5 / (1 - e^2)^1.5; % de_dt = -p / B * n *a * (1 + e) * (1 + 2 * e + e^2)^0.5 / sqrt(1 - e^2); % drp_dt = da_dt * (1 - e) - a * de_dt; dydt = [da_dt; de_dt]; endfunction E = solve_kepler(M, e) % 输入: % M - 平近点角 (弧度) % e - 离心率 (0 ≤ e < 1) % 输出: % E - 偏近点角 (弧度) % 初始猜测 (E0 = M 对小离心率有效) E = M; tol = 1e-12; % 收敛容差 max_iter = 100; % 最大迭代次数 for i = 1:max_iter f = E - e * sin(E) - M; % 方程 f(E) = 0 df = 1 - e * cos(E); % 导数 df/dE delta = f / df; % 牛顿法修正量 E = E - delta; % 更新E % 检查收敛条件 if abs(delta) < tol break; end end % 若未收敛则报错 if i == max_iter warning('未达到收敛容差!'); end end警告: 未达到收敛容差! > 位置:solve_kepler (第 28 行) 位置: dy_dt (第 11 行) 位置: lab1>@(t,y)dy_dt(t,y,B,mu,F_107,Ap) (第 41 行) 位置: ode45 (第 305 行) 位置: lab1 (第 41 行)

给我按文献中的公式中安全系数FOSI的公式(即FOSI=cc1/cc1'=tan(fai1)/tan(fai1')),改变粘聚力(cc1)和内摩擦角(fai1)增加安全系数FOSI,直到使得程序中seigema_anw=0时,此时得到的即为安全系数,求出安全系数的曲线。给出具体思路。程序如下:clc;clear;close all; gamma2=16.1; fai1=deg2rad(38); D=5; C=[0.20:0.25:4]*D; B=pi*D/4; cc1=0.*gamma2*D; Q=0; %地表荷载 delta=0.6.*fai1; %alpha2=deg2rad(71); alpha2=deg2rad(55)+fai1/2; L=D/tan(alpha2); Ka=(1-sin(fai1))./(1+sin(fai1)); %R=D/(2*sqrt(tan(alpha2))); %Ke=1-sin(fai1); %Kp=tan(pi/4+fai1/2).^2; %theta11=atan((Kp-1)+sqrt((Kp-1)^2-4*Kp*tan(delta).^2)./(2*tan(delta)));%Chen %Ka1=tan(pi/4-fai1/2)^2; %Kw=(cos(theta11).^2+Ka1.*sin(theta11).^2)./(1-((1-Ka1).*cos(theta11).^2)/3);%Paik and Salgado(2003)小主应力 % theta2=atan((Kp-1)+sqrt((Kp-1)^2+4*Kp*tan(delta).^2)./(2*tan(delta))) % theta=pi/4+fai1/2; theta=pi/2-1/2.*(asin(sin(delta)./sin(fai1)))+delta./2; theta2=pi/4+fai1/2; beita=0; %Kf=(cos(theta11).^2+Ka1.*sin(theta11).^2)./(((1-Ka1).*sin(theta11).^2)/3+Ka1);%Chen(2015)大主应力 %%%%%楔形体 Kanw1=3.*(cos(theta).^2+Ka.*sin(theta).^2)./(3-(1-Ka).*cos(theta).^2);%23年小主应力 Kanw2=3.*(cos(theta2).^2+Ka.*sin(theta2).^2)./(3-(1-Ka).*cos(theta2).^2); %tao_aw A1=Kanw1.*tan(delta); %tao_s A2=3*(1-Ka)*cos(fai1)/(2*(3-(1-Ka)*cos(theta)^2)); %A2=0; %seigema_s A3=3*(1-Ka)*cos(fai1)^2.*cot(alpha2)/(2*sin(fai1)*(3-(1-Ka)*cos(theta)^2)); %seigema_anw A4=Kanw1.*tan(beita); %tao_aw B1=cc1.*tan(delta).*(1+(1-Ka).*(((Kanw1*cos(theta).^2)/3)-sin(theta).^2))./tan(fai1); %tao_s B2=cc1*(1-Ka)*cos(fai1).^2./2*sin(fai1);%tao_s B2=B2*(1+(1-Ka)*cos(theta)^2./3*(3-(1-Ka)*cos(theta)^2)); %seigema_s B3=cc1*(1-Ka)*cos(fai1).^2.*cot(alpha2)/(2*tan(fai1)*sin(fai1)); B3=B3*(1+(1-Ka)*cos(theta)^2./3*(3-(1-Ka)*cos(theta)^2))-cc1*cot(alpha2)./tan(fai1); %seigema_anw B4=cc1*(1-Ka)*tan(beita).*(((Kanw1*cos(theta).^2)/3)-sin(theta).^2)./tan(fai1); yita_1=A1+A2+A3+A4; yita_2=B1+B2+B3+B4; r=gamma2-2*cc1*sin(alpha2)./B;%对比时把后面的去掉 s=yita_1./(cot(alpha2)+tan(beita)); t=2*sin(alpha2)*Kanw2*tan(fai1)./B; kesi=yita_2./(cot(alpha2)+tan(beita)); %筒仓体 f_11=(Ka-(Ka.*sin(theta2)^2)./3+sin(theta2)^2./3)./(sin(theta2)^2+Ka.*cos(theta2)^2); f_12=f_11.*cos(theta2)^2+sin(theta2)^2./3-1; f_31=(1-Ka).*cos(fai1)./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); f_32=cos(fai1)./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); f_21=(Ka.*(1+sin(fai1))+(1-sin(fai1)))./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); f_22=(cos(theta2)^2-sin(theta2)^2-sin(fai1))./(2.*(sin(theta2)^2+Ka.*cos(theta2)^2)); Kanw0=f_21./f_11; %大主应力侧压力系数(根据东北大学来的) Cpie=L/Kanw0; Ve=pi*B*Cpie*L/4; Se=B*L; C01=min(Ve./Se,C);%筒仓高度 C01_D=C/D; P=2.*cc1.*sqrt(Ka).*(f_32-f_31.*f_12./f_11); E=f_31./f_11; J=2*(B+L)./(B*L); q2=(gamma2-J.*P)./(J.*E)+(Q-(gamma2-J.*P)./(J.*E)).*exp(-J.*E.*C01); figure; hold on; tspan=[0 D-0.01];%楔形体 f=@(z,q) [(r-kesi./(D-z))-((s-1)./(D-z)+t)*q] q0=q2; [z,q]=ode45(f, tspan, q0); seigema_anw=Kanw1*q+cc1*(1-Ka).*(Kanw1*cos(theta).^2/3-sin(theta).^2)./tan(fai1);%水平土压力 S=trapz(z, seigema_anw);%求曲线面积 seigema_t=S./D;%求平均 sei_t=seigema_t./(gamma2.*D); % 绘图 hold on; grid off; box off; set(gca,'LineWidth',1.5,'FontSize',10,'FontName','Times New Roman'); ylabel('$\sigma_{\mathrm{T}} \, /(\mathrm{kPa})$', 'Interpreter', 'latex', 'FontSize', 14, 'FontName', 'Times New Roman'); xlabel('\itC/D', 'FontSize', 14, 'FontName', 'Times New Roman'); ylim([0 0.18]); plot(C01_D, sei_t, '-o', 'Color', 'b', 'LineWidth', 1.5, 'MarkerFaceColor', 'b', 'DisplayName', '本文解'); legend('Location', 'northeast', 'Orientation', 'vertical', 'Interpreter', 'tex', 'FontSize', 12, 'FontName', 'STSong', 'Box', 'off','NumColumns', 2);

最新推荐

recommend-type

二分法和牛顿迭代法求解方程

在上述实验中,二分法求解方程 `1 + x - exp(x) = 0` 的过程中,经过52次迭代得到了满足误差要求的解。 相比之下,牛顿迭代法是一种迭代求解方法,适用于求解可微函数的根。该方法利用函数在某一点的切线与x轴的...
recommend-type

MCP server 项目文件-weather.py

MCP server 项目文件-weather.py
recommend-type

2022版微信自定义密码锁定程序保护隐私

标题《微信锁定程序2022,自定义密码锁》和描述“微信锁定程序2022,自定义密码锁,打开微信需要填写自己设定的密码,才可以查看微信信息和回复信息操作”提及了一个应用程序,该程序为微信用户提供了额外的安全层。以下是对该程序相关的知识点的详细说明: 1. 微信应用程序安全需求 微信作为一种广泛使用的即时通讯工具,其通讯内容涉及大量私人信息,因此用户对其隐私和安全性的需求日益增长。在这样的背景下,出现了第三方应用程序或工具,旨在增强微信的安全性和隐私性,例如我们讨论的“微信锁定程序2022”。 2. “自定义密码锁”功能 “自定义密码锁”是一项特定功能,允许用户通过设定个人密码来增强微信应用程序的安全性。这项功能要求用户在打开微信或尝试查看、回复微信信息时,必须先输入他们设置的密码。这样,即便手机丢失或被盗,未经授权的用户也无法轻易访问微信中的个人信息。 3. 实现自定义密码锁的技术手段 为了实现这种类型的锁定功能,开发人员可能会使用多种技术手段,包括但不限于: - 加密技术:对微信的数据进行加密,确保即使数据被截获,也无法在没有密钥的情况下读取。 - 应用程序层锁定:在软件层面添加一层权限管理,只允许通过验证的用户使用应用程序。 - 操作系统集成:与手机操作系统的安全功能进行集成,利用手机的生物识别技术或复杂的密码保护微信。 - 远程锁定与擦除:提供远程锁定或擦除微信数据的功能,以应对手机丢失或被盗的情况。 4. 微信锁定程序2022的潜在优势 - 增强隐私保护:防止他人未经授权访问微信账户中的对话和媒体文件。 - 防止数据泄露:在手机丢失或被盗的情况下,减少敏感信息泄露的风险。 - 保护未成年人:父母可以为孩子设定密码,控制孩子的微信使用。 - 为商业用途提供安全保障:在商务场合,微信锁定程序可以防止商业机密的泄露。 5. 使用微信锁定程序2022时需注意事项 - 正确的密码管理:用户需要记住设置的密码,并确保密码足够复杂,不易被破解。 - 避免频繁锁定:过于频繁地锁定和解锁可能会降低使用微信的便捷性。 - 兼容性和更新:确保微信锁定程序与当前使用的微信版本兼容,并定期更新以应对安全漏洞。 - 第三方应用风险:使用第三方应用程序可能带来安全风险,用户应从可信来源下载程序并了解其隐私政策。 6. 结语 微信锁定程序2022是一个创新的应用,它提供了附加的安全性措施来保护用户的微信账户。尽管在实施中可能会面临一定的挑战,但它为那些对隐私和安全有更高要求的用户提供了可行的解决方案。在应用此类程序时,用户应谨慎行事,确保其对应用程序的安全性和兼容性有所了解,并采取适当措施保护自己的安全密码。
recommend-type

【自动化脚本提速】:掌握序列生成的5种高效技巧

# 摘要 本文系统地阐述了自动化脚本提速的方法,重点介绍了序列生成的基础理论及其在脚本中的应用。通过探讨不同序列生成方法和高效技巧,本文旨在提高编程效率,优化自动化流程。同时,文中还涉及了高级技术,如嵌套循环、列表推导式和并行处理,这些技术不仅增加了序列生成的复杂性,同时也显著提升了效率。最后,本文通过综合案例分析,展示了一系列序列生成技巧的实际应用,并提出了优化建议和未来研究方向。 #
recommend-type

卷积神经网络中的分层!

<think>我们正在处理一个关于卷积神经网络(CNN)层级结构的问题。用户希望了解CNN的层级结构及其功能。根据提供的引用内容,我们可以整理出以下信息: 1. 引用[1]和[2]指出,一个完整的卷积神经网络通常包括以下层级: - 数据输入层(Input layer) - 卷积计算层(CONV layer) - ReLU激励层(ReLU layer) - 池化层(Pooling layer) - 全连接层(FC layer) - (可能还有)Batch Normalization层 2. 引用[2]详细说明了各层的作用: - 数据输入层:对原始图像
recommend-type

MXNet预训练模型介绍:arcface_r100_v1与retinaface-R50

根据提供的文件信息,我们可以从中提取出关于MXNet深度学习框架、人脸识别技术以及具体预训练模型的知识点。下面将详细说明这些内容。 ### MXNet 深度学习框架 MXNet是一个开源的深度学习框架,由Apache软件基金会支持,它在设计上旨在支持高效、灵活地进行大规模的深度学习。MXNet支持多种编程语言,并且可以部署在不同的设备上,从个人电脑到云服务器集群。它提供高效的多GPU和分布式计算支持,并且具备自动微分机制,允许开发者以声明性的方式表达神经网络模型的定义,并高效地进行训练和推理。 MXNet的一些关键特性包括: 1. **多语言API支持**:MXNet支持Python、Scala、Julia、C++等语言,方便不同背景的开发者使用。 2. **灵活的计算图**:MXNet拥有动态计算图(imperative programming)和静态计算图(symbolic programming)两种编程模型,可以满足不同类型的深度学习任务。 3. **高效的性能**:MXNet优化了底层计算,支持GPU加速,并且在多GPU环境下也进行了性能优化。 4. **自动并行计算**:MXNet可以自动将计算任务分配到CPU和GPU,无需开发者手动介入。 5. **扩展性**:MXNet社区活跃,提供了大量的预训练模型和辅助工具,方便研究人员和开发者在现有工作基础上进行扩展和创新。 ### 人脸识别技术 人脸识别技术是一种基于人的脸部特征信息进行身份识别的生物识别技术,广泛应用于安防、监控、支付验证等领域。该技术通常分为人脸检测(Face Detection)、特征提取(Feature Extraction)和特征匹配(Feature Matching)三个步骤。 1. **人脸检测**:定位出图像中人脸的位置,通常通过深度学习模型实现,如R-CNN、YOLO或SSD等。 2. **特征提取**:从检测到的人脸区域中提取关键的特征信息,这是识别和比较不同人脸的关键步骤。 3. **特征匹配**:将提取的特征与数据库中已有的人脸特征进行比较,得出最相似的人脸特征,从而完成身份验证。 ### 预训练模型 预训练模型是在大量数据上预先训练好的深度学习模型,可以通过迁移学习的方式应用到新的任务上。预训练模型的优点在于可以缩短训练时间,并且在标注数据较少的新任务上也能获得较好的性能。 #### arcface_r100_v1 arcface_r100_v1是一个使用ArcFace损失函数训练的人脸识别模型,基于ResNet-100架构。ArcFace是一种流行的深度学习人脸识别方法,它在损失函数层面上增强类间的区分度。在ArcFace中,通过引入角度余弦的特征分离度,改善了传统的Softmax损失函数,让学习到的人脸特征更加具有鉴别力。 ArcFace的模型文件包括: - model-0000.params: 这是模型权重参数文件。 - model-symbol.json: 这是包含网络结构定义的JSON文件。 #### retinaface-R50 retinaface-R50是基于ResNet-50架构的人脸检测模型,使用RetinaFace框架训练而成。RetinaFace是为了解决传统人脸检测模型在面对小尺寸、遮挡、模糊等复杂情况时识别准确度不高的问题而设计的。它采用一种基于多尺度的金字塔网络结构,能有效处理不同尺度的人脸,并且在特征提取时采用了一种高效的特征融合策略。 Retinaface-R50的模型文件包括: - R50-0000.params: 这是模型权重参数文件。 - R50-symbol.json: 这是包含网络结构定义的JSON文件。 ### 总结 从给定的文件信息中,我们可以看出这些预训练模型是基于MXNet深度学习框架开发的,具有专门针对人脸识别任务的优化。ArcFace模型通过增强特征的区分度,而Retinaface模型通过多尺度处理和高效的特征融合,都展示了在人脸检测和识别方面的先进技术。开发者可以利用这些预训练模型,结合MXNet提供的高级API,快速构建并部署自己的人脸识别解决方案。
recommend-type

【文本处理黑科技】:Shell脚本中序列和数组的高级应用

# 摘要 本文系统地介绍了Shell脚本中序列与数组的基础知识、操作技巧以及在实际应用中的优化方法。第一章和第二章分别对序列和数组的基础知识进行了阐述,重点讲解了序列的生成、操作、文本处理和数组的基本操作与高级特性。第三章探讨了数组在数据处理和内存管理中的性能优势。第四章结合实际案例,展示了序列与数组在日志文件分析、报告生成和配置文件管理
recommend-type

数据库用来干啥的

<think>首先,用户的问题是:“我想了解数据库的用途 请问数据库的主要用途是什么”。这翻译成英文是:“I want to understand the purpose of databases. What are the main uses of databases?” 从系统级指令中,我需要: - 所有行内数学表达式使用 $...$ 格式。 - 独立公式使用 $$...$$ 格式并单独成段。 - LaTeX语法正确。 - 使用中文回答。 - 生成相关问题。 - 回答中引用的段落末尾自然地添加引用标识,例如 [^1]。 - 结构清晰,帮助用户逐步解决问题。 - 回答真实可
recommend-type

微软文字转语音软件发布,支持多国语言高精度识别

微软文字转语音V1.0.zip这个文件包含了一个软件程序,该程序能够实现文字到语音(Text-to-Speech,TTS)的转换,同时也提供了从音频文件到文字的转换功能。下面详细介绍这一技术及其应用的知识点。 首先,文字到语音转换技术,通常称为TTS技术,是计算机科学中的一个重要分支。TTS系统主要由三个关键部分构成:文本分析模块、语音合成模块和声音输出模块。文本分析模块负责理解输入文本的结构和含义,包括分词、词性标注、句法分析、语义分析等。语音合成模块则负责根据文本分析的结果,合成出自然的语音波形。声音输出模块则是将合成的语音波形转换为可听的声音。 微软文字转语音软件具有以下功能和特点: 1. 智能文字和语音互转:该软件能够将用户输入的文字信息转换成自然流畅的语音输出,同时也能将音频文件中的对话转换成文字文本。这种双向转换功能对于多种场景有着重要的应用价值,比如辅助视障人士、提供语音交互界面、制作多语种的语音内容等。 2. 高精度识别各国语言:软件支持高精度的语言识别功能,能处理多种语言的文本或音频。不同语言有不同的语法结构和发音特点,因此支持多语言识别需要对每一种语言都进行深入的研究和算法优化,以确保转换结果的准确性和自然度。 3. 一键拖拽,批量完成:该软件提供简便的操作界面,用户可以通过简单的拖拽动作将需要转换的文本或音频文件直接加入到软件中,进行批量处理。这种操作方式极大地方便了用户,提高了工作效率,尤其在处理大量数据时优势更加明显。 4. 各种音频格式任意选择:用户可以根据需要选择输出的音频格式,比如常见的MP3、WAV等,以便适用于不同的播放设备或播放环境。不同音频格式有其特定的用途,例如MP3格式因为压缩比例高而被广泛用于网络传输和便携式设备,而WAV格式则多用于专业的音频编辑和制作。 软件包中的“resources”文件夹可能包含了支持软件运行的资源文件,如语音合成引擎所需的语音库、语言模型、字典等。而“转换结果”文件夹则可能是软件保存转换后文件的默认位置,用户可以在这里找到转换完成的文字或音频文件。 此外,软件包中的“微软文字转语音V1.0.exe”是一个可执行文件,用户通过运行该文件来启动软件,并使用其提供的各项转换功能。对于IT行业专业人士而言,了解这款软件背后的TTS技术原理和操作逻辑,可以更好地选择合适的解决方案,以满足特定的业务需求。 总结来说,微软文字转语音V1.0.zip中的软件是一款综合性的文字语音转换工具,具有高精度语言识别、高效批量处理、灵活音频格式选择等特点,可以应用于多种场景,提高信息的可访问性和交互性。
recommend-type

【Shell脚本必备】:创建序列的3种方法及高效用法

# 摘要 本文全面探讨了在Shell脚本中创建和优化序列生成的各种方法及其应用场景。首先介绍了序列生成的基本概念和使用基本命令创建序列的技巧,包括for循环、seq命令和算术表达式的应用。随后,深入分析了利用高级Shell特性如数组、复合命令和子shell、以及C风格的for循环来创建复杂序列的技术。文章还探讨了序列在文件批量处理、数据处理分析和自动化脚本中的高效应用。此外,为提升