当前位置：首页 > news >正文

最小作用量原理MATLAB仿真

news 2025/11/10 17:59:09

序

一、自由落体运动的作用量极值推导

1. 定义物理模型

2. 作用量定义与极值条件

3. 代入推导真实轨迹

4. 验证 “作用量最小”

二、MATLAB 仿真：自由落体轨迹的作用量极值可视化

1. 完整仿真代码

2. 仿真结果说明

三、光的折射定律推导（非均匀介质作用量极值）

1. 物理模型定义

2. 作用量表达式与极值条件

3. 代入推导折射定律

四、MATLAB 仿真：光程优化与折射定律可视化

1. 完整仿真代码

2. 仿真结果说明

序

最小作用量原理是物理学核心公理：物理系统的真实运动轨迹，会让 “作用量”（动能与势能的时间积分，记为 S=∫Ldt，L 为拉格朗日量）取极值（通常是最小值）。

它是经典力学、电磁学、相对论、量子力学的统一基础 —— 比如行星轨道、光的折射路径，本质都是 “作用量最小” 的选择，像自然遵循 “最省劲” 的法则。

一、自由落体运动的作用量极值推导

1. 定义物理模型

设质量为 m 的物体从高度 h0 自由下落，重力加速度为 g，忽略空气阻力。

- 动能： $( T = \frac{1}{2}m\dot{y}^2 )$ （$ \dot{y} = \frac{dy}{dt} $ 为竖直速度，向下为正）

- 势能：$ V = -mgy $（以初始位置为势能零点，向下位移 $ y $ 越大，势能越小）

- 拉格朗日量：$ L = T - V = \frac{1}{2}m\dot{y}^2 + mgy $

2. 作用量定义与极值条件

作用量 $ S = \int_{t_0}^{t_1} L dt = \int_{t_0}^{t_1} \left( \frac{1}{2}m\dot{y}^2 + mgy \right) dt $。

根据最小作用量原理，真实运动轨迹 $ y(t) $ 满足**欧拉-拉格朗日方程**：

\[ \frac{d}{dt}\left( \frac{\partial L}{\partial \dot{y}} \right) - \frac{\partial L}{\partial y} = 0 \]

3. 代入推导真实轨迹

- 计算偏导数：$ \frac{\partial L}{\partial \dot{y}} = m\dot{y} $，$ \frac{\partial L}{\partial y} = mg $

- 代入方程：$ \frac{d}{dt}(m\dot{y}) - mg = 0 \implies m\ddot{y} = mg \implies \ddot{y} = g $

- 积分求解：结合初始条件 $ y(0) = 0 $、$ \dot{y}(0) = 0 $，得真实轨迹：

\[ y(t) = \frac{1}{2}gt^2 \]

4. 验证“作用量最小”

假设存在偏离真实轨迹的试探轨迹 $ y'(t) = \frac{1}{2}gt^2 + \delta y(t) $（$ \delta y(t_0) = \delta y(t_1) = 0 $，端点固定），代入作用量公式计算可得 $ S(y'(t)) > S(y(t)) $，证明真实轨迹的作用量取最小值。

二、MATLAB 仿真：自由落体轨迹的作用量极值可视化

1. 完整仿真代码

% 最小作用量原理 - 自由落体轨迹仿真与可视化
clear; clc; close all;%% 1. 物理参数设置
m = 1;          % 物体质量 (kg)
g = 9.8;        % 重力加速度 (m/s²)
t0 = 0;         % 初始时间 (s)
t1 = 2;         % 结束时间 (s)
dt = 0.01;      % 时间步长 (s)
t = t0:dt:t1;   % 时间向量%% 2. 真实轨迹（作用量最小）
y_true = 0.5 * g * t.^2;          % 真实位移 (m)
dot_y_true = g * t;               % 真实速度 (m/s)
L_true = 0.5 * m * dot_y_true.^2 + m * g * y_true;  % 真实拉格朗日量
S_true = trapz(t, L_true);        % 真实作用量（数值积分）%% 3. 试探轨迹（偏离真实轨迹，作用量更大）
% 试探轨迹1：二次函数偏离（幅度0.2m）
y_test1 = 0.5 * g * t.^2 + 0.2 * sin(pi * t / t1);  % 端点固定（sin项在t0/t1为0）
dot_y_test1 = g * t + 0.2 * pi / t1 * cos(pi * t / t1);
L_test1 = 0.5 * m * dot_y_test1.^2 + m * g * y_test1;
S_test1 = trapz(t, L_test1);% 试探轨迹2：线性偏离（幅度0.3m）
y_test2 = 0.5 * g * t.^2 + 0.3 * t .* (t1 - t) / (t1^2/4);  % 中点最大偏离0.3m，端点为0
dot_y_test2 = g * t + 0.3 * (t1 - 2*t) / (t1^2/4);
L_test2 = 0.5 * m * dot_y_test2.^2 + m * g * y_test2;
S_test2 = trapz(t, L_test2);%% 4. 结果可视化
figure('Position', [100, 100, 1200, 800]);% 子图1：位移轨迹对比
subplot(2,2,1);
plot(t, y_true, 'r-', 'LineWidth', 2, 'DisplayName', ['真实轨迹 (S=', num2str(S_true, '%.2f'), ')']);
hold on;
plot(t, y_test1, 'b--', 'LineWidth', 1.5, 'DisplayName', ['试探轨迹1 (S=', num2str(S_test1, '%.2f'), ')']);
plot(t, y_test2, 'g:', 'LineWidth', 1.5, 'DisplayName', ['试探轨迹2 (S=', num2str(S_test2, '%.2f'), ')']);
xlabel('时间 t (s)'); ylabel('位移 y (m)');
title('自由落体轨迹对比（真实轨迹作用量最小）');
legend('Location', 'best'); grid on;% 子图2：拉格朗日量随时间变化
subplot(2,2,2);
plot(t, L_true, 'r-', 'LineWidth', 2, 'DisplayName', '真实轨迹 L');
hold on;
plot(t, L_test1, 'b--', 'LineWidth', 1.5, 'DisplayName', '试探轨迹1 L');
plot(t, L_test2, 'g:', 'LineWidth', 1.5, 'DisplayName', '试探轨迹2 L');
xlabel('时间 t (s)'); ylabel('拉格朗日量 L (J)');
title('拉格朗日量随时间变化');
legend('Location', 'best'); grid on;% 子图3：作用量累积曲线
S_true_cum = cumtrapz(t, L_true);   % 累积作用量
S_test1_cum = cumtrapz(t, L_test1);
S_test2_cum = cumtrapz(t, L_test2);
subplot(2,2,3);
plot(t, S_true_cum, 'r-', 'LineWidth', 2, 'DisplayName', '真实轨迹累积S');
hold on;
plot(t, S_test1_cum, 'b--', 'LineWidth', 1.5, 'DisplayName', '试探轨迹1累积S');
plot(t, S_test2_cum, 'g:', 'LineWidth', 1.5, 'DisplayName', '试探轨迹2累积S');
xlabel('时间 t (s)'); ylabel('累积作用量 S (J·s)');
title('作用量累积曲线（终点为总作用量）');
legend('Location', 'best'); grid on;% 子图4：作用量数值对比柱状图
subplot(2,2,4);
S_values = [S_true, S_test1, S_test2];
labels = {'真实轨迹', '试探轨迹1', '试探轨迹2'};
% 先使用数值索引绘图，再设置标签
bar(1:length(S_values), S_values, 'FaceColor', [0.8, 0.8, 0.8]);
set(gca, 'XTick', 1:length(S_values), 'XTickLabel', labels);  % 手动设置x轴标签
hold on;
text(1:length(S_values), S_values + 0.5, num2str(S_values', '%.2f'), ...'HorizontalAlignment', 'center', 'FontSize', 10);
xlabel('轨迹类型'); ylabel('总作用量 S (J·s)');
title('不同轨迹的总作用量对比');
grid on; ylim([min(S_values)-1, max(S_values)+2]);% 全局标题
sgtitle('最小作用量原理 - 自由落体仿真可视化', 'FontSize', 16, 'FontWeight', 'bold');

2. 仿真结果说明

子图 1：红色实线为真实自由落体轨迹 y(t)=21gt2，蓝色虚线和绿色点线为偏离真实轨迹的试探轨迹，端点均固定在 (0,0) 和 (2,19.6)（t=2s 时真实位移）。
子图 2：拉格朗日量 L=T−V 随时间变化，真实轨迹的 L 变化更平稳。
子图 3：累积作用量曲线，终点值为总作用量，可见真实轨迹的累积作用量始终最小。
子图 4：柱状图直观对比总作用量，真实轨迹的 S 显著小于试探轨迹，验证最小作用量原理。

三、光的折射定律推导（非均匀介质作用量极值）

1. 物理模型定义

设光从折射率为 $ n_1 $ 的介质1（上半空间）入射到折射率为 $ n_2 $ 的介质2（下半空间），分界面为x轴（$ y=0 $）。

- 入射点：$ A(x_1, y_1) $（$ y_1 > 0 $，介质1），出射点：$ B(x_2, y_2) $（$ y_2 < 0 $，介质2）

- 试探入射点：分界面上任意点 $ P(x, 0) $，入射光线 $ AP $ 与法线（y轴）夹角为 $ \theta_1 $，折射光线 $ PB $ 与法线夹角为 $ \theta_2 $

- 作用量（光程）：$ S = n_1 \cdot AP + n_2 \cdot PB $（均匀介质中光程=折射率×几何路程）

2. 作用量表达式与极值条件

- 几何路程：$ AP = \sqrt{(x - x_1)^2 + y_1^2} $，$ PB = \sqrt{(x_2 - x)^2 + y_2^2} $

- 作用量：$ S(x) = n_1\sqrt{(x - x_1)^2 + y_1^2} + n_2\sqrt{(x_2 - x)^2 + y_2^2} $

- 极值条件：作用量取最小值时，对 $ x $ 的导数为0（因只有x一个自由变量，无需欧拉-拉格朗日方程，直接求导）：

\[ \frac{dS}{dx} = 0 \]

3. 代入推导折射定律

- 求导计算：

\[ \frac{dS}{dx} = n_1 \cdot \frac{x - x_1}{\sqrt{(x - x_1)^2 + y_1^2}} + n_2 \cdot \frac{-(x_2 - x)}{\sqrt{(x_2 - x)^2 + y_2^2}} = 0 \]

- 三角函数替换：由几何关系，$ \sin\theta_1 = \frac{x - x_1}{AP} $，$ \sin\theta_2 = \frac{x_2 - x}{PB} $，代入上式得：

\[ n_1 \sin\theta_1 = n_2 \sin\theta_2 \]

即**光的折射定律**，证明真实折射光线满足作用量（光程）最小条件。

四、MATLAB 仿真：光程优化与折射定律可视化

1. 完整仿真代码

% 最小作用量原理 - 光的折射定律仿真（光程优化）
clear; clc; close all;%% 1. 物理参数设置
n1 = 1.0;          % 介质1折射率（空气）
n2 = 1.5;          % 介质2折射率（玻璃）
A = [-2, 1];       % 入射点A (x1, y1)，y1>0（介质1）
B = [2, -1];       % 出射点B (x2, y2)，y2<0（介质2）
x_range = [-3, 3]; % 分界面上试探点P的x范围
x = linspace(x_range(1), x_range(2), 1000); % 试探点x坐标向量%% 2. 计算作用量（光程）与真实入射点
% 作用量（光程）S(x) = n1*AP + n2*PB
AP = sqrt((x - A(1)).^2 + (0 - A(2)).^2);  % 入射路程AP
PB = sqrt((B(1) - x).^2 + (B(2) - 0).^2);  % 折射路程PB
S = n1 * AP + n2 * PB;                    % 总光程（作用量）% 真实入射点：光程最小对应的x（数值求解）
[S_min, idx] = min(S);
x_true = x(idx);  % 真实入射点P的x坐标
P_true = [x_true, 0];  % 真实入射点坐标% 计算入射角θ1和折射角θ2（与法线y轴夹角）
AP_true = sqrt((x_true - A(1))^2 + A(2)^2);
PB_true = sqrt((B(1) - x_true)^2 + B(2)^2);
sin_theta1 = abs(x_true - A(1)) / AP_true;  % 正弦值（绝对值保证角度为正）
sin_theta2 = abs(B(1) - x_true) / PB_true;
theta1 = rad2deg(asin(sin_theta1));         % 入射角（度）
theta2 = rad2deg(asin(sin_theta2));         % 折射角（度）% 验证折射定律：n1*sinθ1 ≈ n2*sinθ2
refraction_law1 = n1 * sin_theta1;
refraction_law2 = n2 * sin_theta2;%% 3. 结果可视化
figure('Position', [100, 100, 1200, 800]);% 子图1：光程（作用量）随试探点x的变化
subplot(2,2,1);
plot(x, S, 'b-', 'LineWidth', 1.5);
hold on;
plot(x_true, S_min, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r', ...'DisplayName', ['最小光程点 (x=', num2str(x_true, '%.2f'), ', S=', num2str(S_min, '%.2f'), ')']);
xlabel('试探入射点x坐标 (m)'); ylabel('光程（作用量）S (m)');
title('光程随分界面入射点的变化（最小值对应真实路径）');
legend('Location', 'best'); grid on;% 子图2：真实折射路径与试探路径对比
subplot(2,2,2);
% 绘制分界面（x轴）
plot(x_range, [0, 0], 'k-', 'LineWidth', 2, 'DisplayName', '介质分界面（y=0）');
hold on;
% 真实路径（AP_true + PB_true）
plot([A(1), P_true(1)], [A(2), P_true(2)], 'r-', 'LineWidth', 2, 'DisplayName', '真实入射光线');
plot([P_true(1), B(1)], [P_true(2), B(2)], 'r-', 'LineWidth', 2);
% 试探路径（以x=0为例）
x_test = 0;
P_test = [x_test, 0];
plot([A(1), P_test(1)], [A(2), P_test(2)], 'g--', 'LineWidth', 1.5, 'DisplayName', '试探入射光线');
plot([P_test(1), B(1)], [P_test(2), B(2)], 'g--', 'LineWidth', 1.5);
% 标记入射点、出射点
scatter(A(1), A(2), 50, 'b', 'filled', 'DisplayName', '入射点A');
scatter(B(1), B(2), 50, 'm', 'filled', 'DisplayName', '出射点B');
scatter(P_true(1), P_true(2), 50, 'r', 'filled', 'DisplayName', '真实入射点P');
scatter(P_test(1), P_test(2), 50, 'g', 'filled', 'DisplayName', '试探入射点P''');
xlabel('x坐标 (m)'); ylabel('y坐标 (m)');
title(['折射定律验证：n1*sinθ1 = ', num2str(refraction_law1, '%.4f'), ', n2*sinθ2 = ', num2str(refraction_law2, '%.4f')]);
axis equal; legend('Location', 'best'); grid on;% 子图3：入射角与折射角可视化
subplot(2,2,3);
% 绘制法线（y轴）
plot([0, 0], [max(A(2), 1.5), min(B(2), -1.5)], 'k:', 'LineWidth', 1.5, 'DisplayName', '法线（y轴）');
hold on;
% 真实路径与角度标注
plot([A(1), P_true(1)], [A(2), P_true(2)], 'r-', 'LineWidth', 2);
plot([P_true(1), B(1)], [P_true(2), B(2)], 'r-', 'LineWidth', 2);
% 角度标注（使用箭头和文本）
annotation('arrow', [0.3, 0.35], [0.7, 0.65], 'Color', 'r');
annotation('textbox', [0.36, 0.68, 0.1, 0.05], 'String', ['θ1 = ', num2str(theta1, '%.1f°')], 'FontSize', 10);
annotation('arrow', [0.6, 0.65], [0.3, 0.25], 'Color', 'r');
annotation('textbox', [0.66, 0.23, 0.1, 0.05], 'String', ['θ2 = ', num2str(theta2, '%.1f°')], 'FontSize', 10);
xlabel('x坐标 (m)'); ylabel('y坐标 (m)');
title('入射角θ1与折射角θ2（相对于法线）');
axis equal; xlim([-1, 1]); ylim([-1.5, 1.5]); grid on;% 子图4：不同折射率组合的折射角变化
subplot(2,2,4);
n1_list = linspace(1.0, 1.2, 20);  % 介质1折射率变化范围
theta2_list = rad2deg(asin(n1_list * sin_theta1 / n2));  % 对应的折射角
plot(n1_list, theta2_list, 'b-', 'LineWidth', 2);
xlabel('介质1折射率n1'); ylabel('折射角θ2 (°)');
title('折射角随介质1折射率的变化（n2固定为1.5）');
grid on;% 全局标题
sgtitle('最小作用量原理 - 光的折射定律仿真可视化', 'FontSize', 16, 'FontWeight', 'bold');