A simulation and analysis of magnetic newton cradle

A multi-variable experiment on magnetic newton cradle and its visualization

function Magenaticpendulumfunc()

% CYPT 2026 Problem 15: 磁性牛顿摆物理仿真定稿版

clear; close all; clc;

%% === 1. 物理参数配置 ===

p.N = 5; % 球数

p.L = 0.45; % 摆长 (m)

p.m = 0.12; % 质量 (kg)

p.d = 0.65; % 悬点间距 (m)

p.g = 9.81;

p.K_mag = 3.2e-1; % 磁耦合强度 (磁偶极子模型)

p.c_damp = 0; % 阻尼系数 F = -cv

%% === 2. 算初始平衡位置 ===

fprintf('正在计算磁力平衡态...\n');

balance_obj = @(th) -p.m * p.g * p.L * sin(th) + calc_mag_torque(th, p);

options_fs = optimoptions('fsolve','Display','off','FunctionTolerance',1e-12);

theta_eq = fsolve(balance_obj, linspace(-0.15, 0.15, p.N)', options_fs);

% 拉开初始球

theta0 = theta_eq;

theta0(1) = theta0(1) - 15 * (pi/180); % 拉开15度释放 %拉n号加theta（n）就好

%theta0(5) = theta0(5) + 15 * (pi/180);

y0 = [theta0; zeros(p.N, 1)];

%% === 3. 高精度动力学求解 ===

t_end = 60; fps = 35;

t_span = linspace(0, t_end, t_end * fps * 8); % 超采样确保平滑

[T, Y] = ode45(@(t,y) physics_kernel(t,y,p), t_span, y0, odeset('RelTol',1e-9));

Theta = Y(:, 1:p.N);

Omega = Y(:, p.N+1:end);

KE = 0.5 * p.m * (p.L * Omega).^2; % 计算系统瞬时动量 [cite: 113]

%% === 5. 物理分析图表 ===

h_anal = figure('Color','w', 'Position', [100 100 1500 900]);

t_lay = tiledlayout(2, 2, 'TileSpacing', 'Compact', 'Padding', 'Normal');

% (A) 动能时空传播图 (窄长比例 + 动能映射)

nexttile;

[X_g, T_g] = meshgrid(1:p.N, T);

pcolor(X_g, T_g, KE);

shading interp; colormap("turbo");

cb = colorbar; cb.Label.String = 'Kinetic Energy E_k (J)';

set(gca, 'XTick', 1:p.N, 'YTick', 0:2:t_end, 'FontSize', 12);

% 关键：强制设置纵横比，使图表窄长清晰

pbaspect([1 2.5 1]);

xlabel('Magnet Ball Index (Discrete)'); ylabel('Physical Time (s)');

title('A. 能量传播路径', 'FontSize', 14);

% (B) 优化步长的相空间轨迹

nexttile;

hold on; colors = lines(p.N);

for n = 1:p.N

plot(Theta(:,n)*(180/pi), Omega(:,n), 'Color', colors(n,:), 'LineWidth', 1.8);

end

grid on; box on;

set(gca, 'XTick', -40:20:40, 'YTick', -3:1:3, 'FontSize', 12);

xlabel('Angle \theta (deg)'); ylabel('Angular Velocity \omega (rad/s)');

title('B. 像空间路径图', 'FontSize', 14);

legend('Ball 1','Ball 2','Ball 3','Ball 4','Ball 5', 'Location', 'northeastoutside');

% (C) 验证扇形平衡态的时间响应图

nexttile([1, 2]);

h_lines = plot(T, Theta * (180/pi), 'LineWidth', 2);

hold on;

for n = 1:p.N

yline(theta_eq(n)*180/pi, '--', 'Color', get(h_lines(n),'Color'), 'Alpha', 0.5, 'HandleVisibility','off');

end

grid on; set(gca, 'XTick', 0:1:t_end, 'YTick', -40:20:40, 'FontSize', 12);

xlabel('Time (s)'); ylabel('Angle (deg)');

title('C. 角速度vs时间', 'FontSize', 14);

fprintf('仿真完成。视频与分析图已保存在当前路径。\n');

%% === 物理子函数 ===

function dydt = physics_kernel(~, y, p)

th = y(1:p.N); om = y(p.N+1:end);

tau_m = calc_mag_torque(th, p);

tau_g = -p.m * p.g * p.L * sin(th);

tau_d = -p.c_damp * (p.L^2) * om;

d_om = (tau_g + tau_m + tau_d) / (p.m * p.L^2);

dydt = [om; d_om];

end

function tau = calc_mag_torque(th, p)

tau = zeros(p.N, 1);

x = (0:p.N-1)' * p.d + p.L * sin(th);

for i = 1:p.N

for j = 1:p.N

if i == j, continue; end

r_ij = x(i) - x(j);

% 磁力建模：1/r^4 偶极子相互作用

f_mag = sign(r_ij) * (p.K_mag / (abs(r_ij)^4 + 1e-10));

tau(i) = tau(i) + f_mag * p.L * cos(th(i));

end

end

end

end