function shangdi5
% 双线摆模型 - 强制左侧线竖直约束（纯无量纲版本）
% ======================== % 1. 无量纲参数定义 % ======================== params = struct(); % 物理参数 params.C0 = 0.4; % 收缩系数 params.tau1_prime = 1; % 光照热弛豫时间（无量纲） params.tau2_prime = 5; % 非光照弛豫时间（无量纲，15/3=5） params.g_prime = 1; % 重力加速度（无量纲，9.8*9/0.3=294） params.c_prime = 100; % 空气阻尼 % 几何参数 params.L_prime = 1/3; % 杆长度（无量纲，0.1/0.3=1/3） params.D_prime = params.L_prime; % A-B距离（无量纲，0.1/0.3=1/3） params.phi_initial = 0; % 初始杆与水平方向夹角（弧度） params.l1_prime = 1; % 左侧线长度（无量纲） params.l_min_prime = 1- params.C0; % LCE最小长度（无量纲，0.3*0.6/0.3=0.6） params.l0_prime = 1; % LCE原长（无量纲） % 其他参数 params.I_prime = 1/12*(params.L_prime)^2; params.alpha =1; % 约束刚度系数 params.beta = 1; % 约束阻尼系数 params.theta0 = pi/500; % 光照角度范围（弧度，约5度） % ======================== % 2. 初始条件设置（无量纲） % ======================== % 关键修改：确保C点在左侧(0, -l1) xc0_prime = params.L_prime / 2; % 质心x = D/2 yc0_prime = -params.l1_prime; % 质心y = -l1 (左侧线竖直) initial_zone = 1; % 初始在光照区 last_switch_time = 0; % 初始切换时间为0 Y0 = [ xc0_prime; % 质心x（无量纲） yc0_prime; % 质心y（无量纲） params.phi_initial; % 初始角度 0; 0; 0; % 初始速度 last_switch_time; % 最后切换时间 initial_zone; % 当前区域（1=光照，0=非光照） 1; % 切换时的LCE长度（初始为1） ]; % ======================== % 3. 数值求解 % ======================== tspan_prime = [0, 500]; % 无量纲时间范围 dt_prime = 0.01; % 无量纲时间步长 num_steps = floor((tspan_prime(2)-tspan_prime(1))/dt_prime) + 1; t_prime = linspace(tspan_prime(1), tspan_prime(2), num_steps)'; % 初始化状态数组 Y = zeros(num_steps, length(Y0)); Y(1,:) = Y0'; % 初始化验证：打印初始C点坐标 xc_val = Y(1,1); yc_val = Y(1,2); phi_val = Y(1,3); A_x = xc_val - params.L_prime/2 * cos(phi_val); A_y = yc_val - params.L_prime/2 * sin(phi_val); fprintf('初始C点坐标: (%.4f, %.4f)\n', A_x, A_y); % 龙格-库塔法求解（带区域切换检测） for i = 1:num_steps-1 current_time = t_prime(i); current_state = Y(i,:)'; % 在调用ODE函数前检测区域切换 [new_state, switched] = check_zone_switch(current_state, current_time, params); if switched Y(i,:) = new_state'; % 更新当前状态 current_state = new_state; % 使用更新后的状态 end % 龙格-库塔迭代 k1 = double_pendulum_ode_non_dim(current_time, current_state, params); k2 = double_pendulum_ode_non_dim(current_time + dt_prime/2, current_state + (dt_prime/2)*k1, params); k3 = double_pendulum_ode_non_dim(current_time + dt_prime/2, current_state + (dt_prime/2)*k2, params); k4 = double_pendulum_ode_non_dim(current_time + dt_prime, current_state + dt_prime*k3, params); next_state = current_state + (dt_prime/6)*(k1 + 2*k2 + 2*k3 + k4); Y(i+1,:) = next_state'; end % ======================== % 4. 结果提取 % ======================== phi = Y(:,3); % 杆与水平方向夹角 dphi = Y(:,6); % 角速度 xc_prime = Y(:,1); % 质心x坐标 yc_prime = Y(:,2); % 质心y坐标 % 计算LCE长度 l2_prime = zeros(size(t_prime)); for i = 1:length(t_prime) [l2_prime(i), ~, ~] = compute_lce_length(t_prime(i), Y(i,:)', params); end % LCE与竖直方向夹角计算 theta2 = zeros(size(phi)); for i = 1:length(t_prime) xc_val = Y(i,1); yc_val = Y(i,2); phi_val = Y(i,3); B_x = xc_val + params.L_prime/2*cos(phi_val); B_y = yc_val + params.L_prime/2*sin(phi_val); dx = B_x - params.D_prime; dy = B_y; theta2(i) = atan2(dx, -dy); end theta2_deg = theta2 * 180/pi; % 计算普通线与竖直方向夹角θ1 theta1 = zeros(size(phi)); for i = 1:length(t_prime) xc_val = Y(i,1); yc_val = Y(i,2); phi_val = Y(i,3); % C点坐标（左侧端点） C_x = xc_val - params.L_prime/2*cos(phi_val); C_y = yc_val - params.L_prime/2*sin(phi_val); % 计算普通线与竖直方向的夹角（逆时针为正） theta1(i) = atan2(C_x, -C_y); end theta1_deg = theta1 * 180/pi; % ======================== % 5. 图形显示（无量纲） % ======================== % 图1：杆的倾斜角随时间变化 figure; plot(t_prime, phi, 'b-', 'LineWidth', 1.2); xlabel('无量纲时间'); ylabel('\phi (rad)'); title('杆与水平方向夹角随时间变化'); grid on; % 图2：LCE纤维长度变化 figure; plot(t_prime, l2_prime, 'r-', 'LineWidth', 1.2); hold on; plot(t_prime, ones(size(t_prime)), 'k--', 'LineWidth', 1); % 原长参考线 hold off; xlabel('无量纲时间'); ylabel('无量纲LCE纤维长度'); title('LCE纤维长度随时间变化'); grid on; % 图3：角度与角速度的相图 figure; time_mask = t_prime >= 400 & t_prime <= 500; plot(phi(time_mask), dphi(time_mask), 'g-', 'LineWidth', 1.2); xlabel('\phi (rad)'); ylabel('无量纲角速度 (rad/\tau_1)'); title('400-500无量纲时间内角度与角速度的相图'); grid on; % 图4：LCE与竖直方向夹角 figure; plot(t_prime, theta2, 'm-', 'LineWidth', 1.2); hold on; plot(t_prime, ones(size(t_prime))*params.theta0, 'k--'); plot(t_prime, -ones(size(t_prime))*params.theta0, 'k--'); hold off; xlabel('无量纲时间'); ylabel('\theta_2 (rad)'); title('LCE与竖直方向夹角随时间变化'); grid on; % 图5：普通线与竖直方向夹角 figure; plot(t_prime, theta1, 'c-', 'LineWidth', 1.2); xlabel('无量纲时间'); ylabel('\theta_1 (rad)'); title('普通线与竖直方向夹角随时间变化'); grid on; % 图6：质心坐标随时间变化 figure; subplot(2,1,1); plot(t_prime, xc_prime, 'b-', 'LineWidth', 1.2); ylabel('x_c'); title('质心x坐标随时间变化'); grid on; subplot(2,1,2); plot(t_prime, yc_prime, 'r-', 'LineWidth', 1.2); xlabel('无量纲时间'); ylabel('y_c'); title('质心y坐标随时间变化'); grid on; % 图7：运动动画 figure; ax = gca; axis equal; xlim([-0.5, max(xc_prime)+0.5]); ylim([min(yc_prime)-0.5, 0.5]); xlabel('无量纲X'); ylabel('无量纲Y'); title('双线摆运动动画'); grid on; hold on; % 绘制横梁与光照区域 plot([-0.5, params.D_prime + 0.5], [0, 0], 'k-', 'LineWidth', 3); % 横梁 % 光照区域（以B点为中心） light_radius = 0.2; theta_range = linspace(-params.theta0, params.theta0, 100); light_x = params.D_prime + light_radius * sin(theta_range); light_y = -light_radius * cos(theta_range); fill(light_x, light_y, 'y', 'FaceAlpha', 0.4, 'EdgeColor', 'none'); % 初始化动画元素 hLine1 = plot([0,0], [0,0], 'b-', 'LineWidth', 2); % 左侧线（蓝色） hLine2 = plot([0,0], [0,0], 'r-', 'LineWidth', 2); % LCE线（红色） hRod = plot([0,0], [0,0], 'k-', 'LineWidth', 4); % 杆 hPointA = plot(0, 0, 'bo', 'MarkerSize', 8, 'MarkerFaceColor', 'b'); % A点 hPointB = plot(0, 0, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r'); % B点 hCenter = plot(0, 0, 'go', 'MarkerSize', 8, 'MarkerFaceColor', 'g'); % 质心 hLightText = text(params.D_prime - 0.3, -0.3, '光照中', 'Color', 'red', 'FontSize', 14, 'Visible', 'off'); hold off; % 动画循环 for i = 1:5:num_steps if ~ishandle(gcf), break; end xc_val = Y(i,1); yc_val = Y(i,2); phi_val = Y(i,3); % C点（左侧端点） A_x = xc_val - params.L_prime/2*cos(phi_val); A_y = yc_val - params.L_prime/2*sin(phi_val); % D点（右侧端点） B_x = xc_val + params.L_prime/2*cos(phi_val); B_y = yc_val + params.L_prime/2*sin(phi_val); % 更新图形 set(hLine1, 'XData', [0, A_x], 'YData', [0, A_y]); % 左侧线 set(hLine2, 'XData', [params.D_prime, B_x], 'YData', [0, B_y]); % LCE线 set(hRod, 'XData', [A_x, B_x], 'YData', [A_y, B_y]); % 杆 set(hPointA, 'XData', A_x, 'YData', A_y); % C点 set(hPointB, 'XData', B_x, 'YData', B_y); % D点 set(hCenter, 'XData', xc_val, 'YData', yc_val); % 质心 % 光照状态显示 current_theta = atan2(B_x - params.D_prime, -B_y); if current_theta >= -params.theta0 && current_theta <= params.theta0 set(hLightText, 'Visible', 'on'); else set(hLightText, 'Visible', 'off'); end drawnow; pause(0.01); end end function [new_state, switched] = check_zone_switch(current_state, t, params) % 检测区域切换并更新状态 new_state = current_state; switched = false; % 解包当前状态 xc = current_state(1); yc = current_state(2); phi = current_state(3); current_zone = current_state(8); % 计算当前角度theta2 B_x = xc + params.L_prime/2*cos(phi); B_y = yc + params.L_prime/2*sin(phi); dx = B_x - params.D_prime; dy = B_y; theta2 = atan2(dx, -dy); % 弧度 % 判断当前应属区域 if theta2 >= -params.theta0 && theta2 <= params.theta0 new_zone = 1; % 光照区 else new_zone = 0; % 暗区 end % 检查是否切换 if new_zone ~= current_zone % 计算切换时刻的LCE长度 [l_switch, ~, ~] = compute_lce_length(t, current_state, params); % 更新状态 new_state(7) = t; % 更新最后切换时间 new_state(8) = new_zone; % 更新区域标志 new_state(9) = l_switch; % 更新切换时刻的LCE长度 switched = true; end end function [l2_prime, dl2dt, d2l2dt2] = compute_lce_length(t, Y, params) % 计算LCE长度及其导数 last_switch_time = Y(7); current_zone = Y(8); l_switch = Y(9); % 切换时刻的长度 local_time = t - last_switch_time; % 当前区域的局部时间 if current_zone == 1 % 光照区 % 光照区目标长度：最小长度 l_target = params.l_min_prime; l2_prime = l_switch + (l_target - l_switch) * (1 - exp(-local_time)); dl2dt = (l_target - l_switch) * exp(-local_time); d2l2dt2 = -(l_target - l_switch) * exp(-local_time); else % 暗区 % 暗区目标长度：原始长度 l_target = 1; l2_prime = l_switch + (l_target - l_switch) * (1 - exp(-local_time / params.tau2_prime)); dl2dt = (l_target - l_switch) * (1/params.tau2_prime) * exp(-local_time / params.tau2_prime); d2l2dt2 = - (l_target - l_switch) * (1/(params.tau2_prime^2)) * exp(-local_time / params.tau2_prime); end end function dYdt = double_pendulum_ode_non_dim(t, Y, params) % 无量纲化的ODE求解函数（加强约束版） Y = Y(:); % 解包参数 L_prime = params.L_prime; g_prime = params.g_prime; D_prime = params.D_prime; l1_prime = params.l1_prime; c_prime = params.c_prime; I_prime = params.I_prime; alpha = params.alpha; beta = params.beta; % 计算LCE长度及其导数 [l2_prime, dl2dt, d2l2dt2] = compute_lce_length(t, Y, params); % 解包状态变量 xc = Y(1); yc = Y(2); phi = Y(3); dxc = Y(4); dyc = Y(5); dphi = Y(6); % 端点坐标 A_x = xc - L_prime/2*cos(phi); A_y = yc - L_prime/2*sin(phi); B_x = xc + L_prime/2*cos(phi); B_y = yc + L_prime/2*sin(phi); % 关键约束1：左侧线长度恒定 f1 = A_x^2 + A_y^2 - l1_prime^2; % 关键约束2：左侧线竖直（A点x坐标=0） f2 = A_x; % 新增强制约束 % LCE线约束 f3 = (B_x - D_prime)^2 + B_y^2 - l2_prime^2; % 导数计算 dA_x = dxc + (L_prime/2)*sin(phi)*dphi; dA_y = dyc - (L_prime/2)*cos(phi)*dphi; dB_x = dxc - (L_prime/2)*sin(phi)*dphi; dB_y = dyc + (L_prime/2)*cos(phi)*dphi; df1dt = 2*A_x*dA_x + 2*A_y*dA_y; df2dt = dA_x; % f2的导数 df3dt = 2*(B_x - D_prime)*dB_x + 2*B_y*dB_y - 2*l2_prime*dl2dt; % 雅可比项 term1 = (L_prime/2) * (A_x*sin(phi) - A_y*cos(phi)); term2 = (L_prime/2) * (-(B_x - D_prime)*sin(phi) + B_y*cos(phi)); % 右端项 R1 = -c_prime * dxc; % 水平阻尼力 R2 = -g_prime - c_prime * dyc; % 竖直方向力 R3 = -1/4*c_prime * (L_prime)^2*dphi; % 旋转阻尼力矩 % 约束方程的右端项 R4 = -2*(dA_x^2 + dA_y^2) - alpha*df1dt - beta^2*f1 ... - L_prime*dphi^2*(A_x*cos(phi) + A_y*sin(phi)); R5 = -2*(dA_x^2) - alpha*df2dt - beta^2*f2; % f2约束的右端项 R6 = -2*(dB_x^2 + dB_y^2) - alpha*df3dt - beta^2*f3 ... + L_prime*dphi^2*((B_x - D_prime)*cos(phi) + B_y*sin(phi)) ... + 2*l2_prime*d2l2dt2 + 2*dl2dt^2; % 构建矩阵方程（6x6系统） A_mat = zeros(6,6); A_mat(1:3,4) = [1; 0; 0]; A_mat(1:3,5) = [0; 1; 0]; A_mat(1:3,6) = [0; 0; I_prime]; % 约束梯度项 A_mat(1:3,1) = [ -2*A_x; -2*A_y; -2*term1 ]; A_mat(1:3,2) = [ -2*(B_x-D_prime); -2*B_y; -2*term2 ]; % f2约束的雅可比 A_mat(4:6,1) = [2*A_x; 1; 2*(B_x-D_prime)]; A_mat(4:6,2) = [2*A_y; 0; 2*B_y]; A_mat(4:6,3) = [2*term1; 0; 2*term2]; B_vec = [R1; R2; R3; R4; R5; R6]; % 求解线性系统 sol = A_mat \ B_vec; % 状态导数向量 dYdt = zeros(9,1); dYdt(1) = Y(4); % dxc/dt = v_xc dYdt(2) = Y(5); % dyc/dt = v_yc dYdt(3) = Y(6); % dphi/dt = omega dYdt(4) = sol(1); % d(v_xc)/dt dYdt(5) = sol(2); % d(v_yc)/dt dYdt(6) = sol(3); % d(omega)/dt dYdt(7) = 0; % d(last_switch_time)/dt = 0 dYdt(8) = 0; % d(current_zone)/dt = 0 dYdt(9) = 0; % d(l_switch)/dt = 0 end 这是一篇matlab代码，代码中用了加速度级约束和惩罚函数，导致改变阻尼系数对结果不产生任何影响。请用牛顿迭代法来处理约束，并计算约束乘子，请给出完整·代码
最新发布
扫一扫
300
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
