Dimensional Analysis量纲分析入门

最新推荐文章于 2026-06-18 14:48:02 发布

原创最新推荐文章于 2026-06-18 14:48:02 发布 · 1.1k 阅读

23 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

标签

#matlab #dynamics #dimensional #量纲分析 #动力学

马特拉博Matlab 专栏收录该内容

51 篇文章

订阅专栏

博客文章：Dimensional Analysis量纲分析入门

幂律（Power Law）

在物理学中，幂律（Power Law）是指量之间的幂次关系，通常表示为：
$kx_1^{n_1}x_2^{n_2}\cdots x_m^{n_m}$

其中， $k$ 是常数， $x_1, x_2, \cdots, x_m$ 是变量， $n_1, n_2, \cdots, n_m$ 是幂次。

定理：有单位的物理量，只能形成幂律关系。

这个定理的证明很简单（！）。

在这里插入图片描述

证明：不失一般性，假设 $y$ 和 $x$ 是两个有单位的物理量，那么 $y$ 和 $x$ 的关系可以表示为：

$y = f (x)$

量纲，也就是物理单位，通常本身也是一种比例关系，例如长度，很原始的基准是手肘到中指尖的长度（英尺），很科幻的现代基准是光在真空中1/299792458秒内走过的距离（米）。有物理单位的量，可以通过取比值的方式来消去单位。

$\frac{y_1}{y_2} = \frac{f(x_1)}{f(x_2)}$

这个式子，就是一个与单位无关的公式，也就是说，无论 $x$ 的单位是什么， $y$ 的单位是什么，这个式子都成立。

那么这里就有一个决定性的推理。

$\frac{f(x_1)}{f(x_2)} = \frac{\alpha x_1}{\alpha x_2}$

也就是， $x$ 取不同的单位，体现为单位转换系数 $\alpha$ ， $\frac{y_1}{y_2}$ 不变。只要这个式子成立，很容易就能得到： $y$ 和 $x$ 就满足幂律关系。

对上面的式子求取 $\alpha$ 的导数，可以得到：

$x_1f(\alpha x_2) f'(\alpha x_1) = x_2 f(\alpha x_1) f'(\alpha x_2)$

对所有的 $\alpha, x_1, x_2$ ，这个式子都成立，我们取 $\alpha = 1, x_2 = 1, x_1 = x$ ，可以得到：

$\frac{x f'(x)}{f(x)} = \frac{f'(1)}{f(1)} \equiv k$

这里的 $k$ 是一个常数。

$\int \frac{f'(x)}{f(x)} = \int \frac{k}{x}$

计算不定积分，可以得到：

$\ln f(x) = k \ln x + c$

取指数，可以得到：

$f(x) = C x^k$

这里 $C$ 是另外一个积分常数。Quod Erat Demonstrandum （Q.E.D.）。

当然，从一般的物理规律出发，也能理解为什么有单位的量只能出现在幂律关系中。就比如一个长度量 $l$ ， $e^l$ 的含义是什么？ $\sin(l)$ 的含义是什么？从数学上来看

$e^l = 1 + l + \frac{l^2}{2!} + \frac{l^3}{3!} + \cdots$

等于把无量纲量、长度量、面积量、体积量……全部加在一起，这显然是没有任何物理意义的。

SI量纲

SI量纲，也就是国际单位制，是国际上通用的物理量单位体系。SI量纲有7个基本量，参考国标GB 3100-1993的3.1节SI基本单位：

长度（ $L$ ）, 单位：米（ $m$ ）
质量（ $M$ ）, 单位：千克（ $k g$ ）
时间（ $T$ ）, 单位：秒（ $s$ ）
电流（ $I$ ）, 单位：安培（ $A$ ）
热力学温度（ $Θ$ ）, 单位：开尔文（ $K$ ）
物质的量（ $N$ ）, 单位：摩尔（ $m o l$ ）
发光强度（ $J$ ）, 单位：坎德拉（ $c d$ ）

我们把一个物理量的单位信息记为 $[y]$ ，表达为上述基本量的幂次形式，例如速度的单位是 $\frac{m}{s}$ ，那么 $v] = L T^{-1}$ ，加速度的单位是 $\frac{m}{s^2}$ ，那么 $a] = L T^{-2}$ 。

所以对于有单位物理量 $y$ ，可以表示为：

$a^\alpha b^\beta c^\gamma \cdots$

其中 $\cdots$ 是有单位的物理参数， $C$ 是常数， $\alpha, \beta, \gamma, \cdots$ 是幂次。

$[a]^\alpha [b]^\beta [c]^\gamma \cdots$

再结合上面的SI量纲，可以得到一组方程，分别对应7个基本量：

$\begin{split} &[y] = L^{\alpha_y} M^{\beta_y} T^{\gamma_y} I^{\delta_y} Θ^{\epsilon_y} N^{\zeta_y} J^{\eta_y} \\ &[a] = L^{\alpha_a} M^{\beta_a} T^{\gamma_a} I^{\delta_a} Θ^{\epsilon_a} N^{\zeta_a} J^{\eta_a} \\ &[b] = L^{\alpha_b} M^{\beta_b} T^{\gamma_b} I^{\delta_b} Θ^{\epsilon_b} N^{\zeta_b} J^{\eta_b} \\ &[c] = L^{\alpha_c} M^{\beta_c} T^{\gamma_c} I^{\delta_c} Θ^{\epsilon_c} N^{\zeta_c} J^{\eta_c} \\ &\cdots \end{split}$

这样可以得到一组方程：

$\begin{split} &\alpha_y = \alpha_a \cdot \alpha + \alpha_b \cdot \beta + \alpha_c \cdot \gamma + \cdots \\ &\beta_y = \beta_a \cdot \alpha + \beta_b \cdot \beta + \beta_c \cdot \gamma + \cdots \\ &\gamma_y = \gamma_a \cdot \alpha + \gamma_b \cdot \beta + \gamma_c \cdot \gamma + \cdots \\ &\delta_y = \delta_a \cdot \alpha + \delta_b \cdot \beta + \delta_c \cdot \gamma + \cdots \\ &\epsilon_y = \epsilon_a \cdot \alpha + \epsilon_b \cdot \beta + \epsilon_c \cdot \gamma + \cdots \\ &\zeta_y = \zeta_a \cdot \alpha + \zeta_b \cdot \beta + \zeta_c \cdot \gamma + \cdots \\ &\eta_y = \eta_a \cdot \alpha + \eta_b \cdot \beta + \eta_c \cdot \gamma + \cdots \\ \end{split}$

举个例子

我们来分析一个单摆的周期。单摆的长度为 $l$ ，质量为 $m$ ，角频率为 $\omega$ ，重力加速度为 $g$ 。

在这里插入图片描述

我们写成：

$\omega = C l^\alpha m^\beta g^\gamma$

量纲分析有：

$\begin{split} &[l] = L \\ &[m] = M \\ &[g] = L T^{-2} \\ &[\omega] = T^{-1} \end{split}$

因此有：

$T^{-1} = L^\alpha M^\beta (L T^{-2})^\gamma$

$\left\{ \begin{split} &\alpha + \gamma = 0 \\ &\beta = 0 \\ & -2\gamma = -1 \end{split} \right.$

解得：

$\left\{ \begin{split} &\alpha = -\frac{1}{2} \\ &\beta = 0 \\ &\gamma = \frac{1}{2} \end{split} \right.$

所以有：

$\omega = C l^{-\frac{1}{2}} g^{\frac{1}{2}} = C \sqrt{\frac{g}{l}}$

实际上，通过动力学积分，也能够很简单地得到：

$\omega = \sqrt{\frac{g}{l}} = 2\pi f$

频率为
$\frac{1}{2\pi} \sqrt{\frac{g}{l}}$

周期为

$\frac{1}{f} = 2\pi \sqrt{\frac{l}{g}}$

这个很简单的推导就留给读者作为练习（！）。

在这里插入图片描述

进一步的讨论

当然，前面的7个方程，其解的情况可以分为三类：

无解：这个是时候，我们就知道肯定丢失了重要的参数
唯一解：非常幸运，只需要通过实验确定常数 $C$
多解：依然得到了一些有用的信息，可以用于指导实验

对于无解和唯一解的情形，我们很容易理解。那么对于多解的情形，物理上的含义如何呢？我们如何选择呢？有什么意义呢？

我们再考虑一个例子，用速度 $V$ 把一个质量为 $m$ 的物体向上抛出，这个球所受到的空气阻力服从平方率，也就是阻力 $k V^2$ 。这里的 $k$ 是阻力系数。问，这个物体上升的高度 $h$ 和什么有关？

$m^\alpha g^\beta k^\gamma V^\delta$

有 $k] = [force]/[velocity]^2 = M L T^{-2} L^{-2} / (L T^{-1})^2 = M L^{-1}$ ，所以有：

$\left\{ \begin{split} & \alpha + \gamma = 0 \\ & \beta - \gamma + \delta = 1\\ & -2\beta -\delta = 0 \end{split} \right.$

这个方程组有无穷多解。例如前面所说的

$\left\{ \begin{split} & \alpha = 1 \\ & \beta = 0 \\ & \gamma = -1 \\ & \delta = 0 \end{split} \right.$

这就表明，上升高度与 $m / k$ 成正比，如果我们把 $h$ 替换为 $h / (m / k)$ ，就可以得到一个无量纲量。进行同样的量纲分析

$\frac{h}{m/k} = C m^\alpha g^\beta k^\gamma V^\delta$

得到如下的方程，因为是 $h / (m / k)$ 是无量纲量，所以方程的右边全部是0。

$\left\{ \begin{split} & \alpha + \gamma = 0 \\ & \beta - \gamma + \delta = 0\\ & -2\beta -\delta = 0 \end{split} \right.$

解得：

$\left\{ \begin{split} & \beta = \alpha \\ & \gamma = -\alpha \\ & \delta = -2\alpha \end{split} \right.$

所以有：

$\frac{h}{m/k} = C m^\alpha g^\alpha k^{-\alpha} V^{-2\alpha} = C \left(\frac{mg}{k V^2}\right)^\alpha$

从这里我们可以得出结论， $\lambda = \frac{mg}{k V^2}$ 是我们问题中的唯一需要考虑的无量纲量。

$\frac{h}{m/k} = f \left(\frac{mg}{k V^2}\right) \text{, i.e. } h = \frac{m}{k} f \left(\frac{mg}{k V^2}\right)$

因为 $\lambda$ 是一个无量纲量，所以这个函数 $f$ 就不见得是幂律关系。实际上，如果我们简单（！）地分析动力学，可以得到：

$f(\lambda) = \frac{1}{2} \ln (1+\lambda ^{-1})$

推导到这里就有一个问题，因为前面我们选择的组合是随意的，所以的，这个结果有意义吗？这里只是简单展示一下。

例如，

$\left\{ \begin{split} & \alpha + \gamma = 0 \\ & \beta - \gamma + \delta = 1\\ & -2\beta -\delta = 0 \end{split} \right.$

的解还可能是：

$\left\{ \begin{split} & \alpha = 0 \\ & \beta = -1 \\ & \gamma = 0 \\ & \delta = 2 \end{split} \right.$

这个时候，上面的过程就变成了：

$\frac{h}{V^2/g} = C m^\alpha g^\beta k^\gamma V^\delta$

当然，这次的量纲方程组求解会得到同样的结果，也就是 $\lambda = \frac{mg}{k V^2}$ 。

$\frac{V^2}{g} \hat{f} \left(\frac{mg}{k V^2}\right)$

观察发现，

$\tilde{f}(\lambda) = \hat{f} / \lambda$

$\frac{m}{k} \tilde{f} \left(\frac{mg}{k V^2}\right)$

因此，选择不同的解，得到的结果是等价的。

通常，我们也把这个公式写成：

$g\left(\frac{mg}{kV^2}, \frac{kh}{m}\right) = 0$

这里 $g$ 是一个函数， $\frac{mg}{kV^2}$ 和 $\frac{kh}{m}$ 是两个无量纲量。当然， $g$ 可以有很多种用 $f$ 表示的形式。

$g (x, y) = y - f (x)$

$g(x, y) = y^2 - f^2(x)$

或者

$\ln \frac{1+y}{1+f(x)}$

如果在上面的例子里面，我们再增加变量，就是丢球的角度 $\theta$ ，因为 $\theta$ 是角度，一个无量纲量，那么我们就能写成：

$\frac{m}{k} F \left(\frac{mg}{k V^2}, \theta\right) = \frac{m}{k} F \left(\lambda, \theta\right)$

实际上， $F$ 无法写出解析形式。

那么，我们这么做有什么意义呢？

通过量纲分析，我们把一个依赖于多个物理量的函数，简化为一个依赖于较少的无量纲量的函数，这非常有助于我们对问题的理解，并且可以指导我们进行实验设计。

抛石问题分析

考虑一个二维的抛石问题，定义位置向量为 $\vec{x} = (x, y)$ ，抛石角度为 $\theta$ ，抛石初速度为 $V$ 。

$\ddot{\vec{x}} = m\vec{g} - k |\dot{\vec{x}}| \dot{\vec{x}}$

展开得到：

$\begin{split} & m \ddot{x} = -k \sqrt{\dot{x}^2 + \dot{y}^2} \dot{x} \\ & m \ddot{y} = -k \sqrt{\dot{x}^2 + \dot{y}^2} \dot{y} - mg \end{split}$

$t = 0$ 时， $\vec{x} = (0, 0)$ ， $\dot{\vec{x}} = (V \cos \theta, V \sin \theta)$ 。

经过上面的量纲分析，我们可以发现，长度的量纲组合是 $m / k$ ，速度的量纲组合是 $V$ ，那么时间量纲组合是 $m / kV$ 。我们就可以按照量纲组合的把量纲量转化成无量纲量。

$\frac{x}{m/k}, \quad Y = \frac{y}{m/k}, \quad T = \frac{t}{m/kV}$

根据链式法则，

$\dot{x} = \frac{dT}{dt}\frac{dx}{dT} = \frac{kV}{m} \frac{d}{dT}\left(\frac{m}{k} X\right) = V X'$

$\ddot{x} = \frac{dT}{dt}\frac{d\dot{x}}{dT} = \frac{kV}{m} \frac{d}{dT}(V X') = \frac{kV^2}{m} X''$

最终可以得到：

$\begin{split} & kV^2 X'' = -k \sqrt{V^2 X'^2 + V^2 Y'^2} V X' \\ & kV^2 Y'' = -k \sqrt{V^2 X'^2 + V^2 Y'^2} V Y' - mg \end{split}$

根据 $\lambda = \frac{mg}{kV^2}$ ，可以得到：

$\begin{split} & X'' = - \sqrt{X'^2 + Y'^2} X' \\ & Y'' = - \sqrt{X'^2 + Y'^2} Y' - \lambda \end{split}$

$T = 0$ 时， $\cos \theta, Y' = \sin \theta$ 。这个方程组的初始值和方程都依赖于无量纲量，整个问题只依赖于两个参数 $\lambda$ 和 $\theta$ 。对于其运动的最高点（对应时间 $T_0$ ）， $Y'(T_0) = 0$ ，此时 $h = (m/k) Y(T_0)$ 。

这就成功地把一个 $\theta)$ 的问题，转化为了一个 $(\lambda, \theta)$ 的求解。我们这里甚至可以给 $\lambda$ 一个名称，考虑到这个参数实际上描述了重力与阻力的相对关系，我们就叫它格拉维-夸德雷特-埃尔雷兹斯唐斯数，简称Grr数。

程序实现

都已经写到这里了，很难不搞一个程序。

function [value, isterminal, direction] = heightEvent(t, y, Y0)
% 当Y回到初始高度时停止
value = y(2) - Y0;      % 当前高度与初始高度的差
isterminal = 1;         % 停止积分
direction = -1;         % 只在下降穿过初始高度时停止
end

function dydt = dimensionlessODE(t, y, lambda)
% y = [X; Y; X'; Y']
X = y(1);
Y = y(2);
Xp = y(3);
Yp = y(4);

% Calculate velocity magnitude
V = sqrt(Xp^2 + Yp^2);

% ODE system
dydt = zeros(4,1);
dydt(1) = Xp;
dydt(2) = Yp;
dydt(3) = -V * Xp;
dydt(4) = -V * Yp - lambda;
end

这两个函数定义了ODE和终止条件（高度再次为0的事件）。

调用ODE45求解，得到结果。

        function [T, X, Y] = solveDimensionlessODE(app)
            % 初始条件
            X0 = 0;
            Y0 = 0;
            Xp0 = cos(app.Theta);
            Yp0 = sin(app.Theta);
            
            % 求解ODE直到Y回到初始高度
            y0 = [X0; Y0; Xp0; Yp0];
            options = odeset('Events', @(t,y) heightEvent(t,y,Y0), ...
                'RelTol', 1e-8, ...  % 相对误差
                'AbsTol', 1e-8, ...  % 绝对误差
                'MaxStep', 0.01);    % 最大步长
            [T, Y] = ode45(@(t,y) dimensionlessODE(t, y, app.Lambda), [0 100], y0, options);
            
            X = Y(:,1);
            Y = Y(:,2);
        end

界面什么的就不纠结了，左边是几个物理参数：质量、阻力系数、重力加速度、初始速度、抛射角度；下面是相应的参考长度、参考时间、最大投射高度（有量纲、无量纲）。

app = ProjectileApp;
exportapp(app.Figure, "mainui.png")

在这里插入图片描述

右边两个标签，分别是有量纲的轨迹和两个无量纲参数下最大射程的图形。

我们就演示一个东西，就是不同的质量下的投射距离。我们一通调节质量，得到如下的结果：

exportgraphics(app.NonDimAxes, "ndt.png")

在这里插入图片描述

exportgraphics(app.DimensionalAxes, "dt.png")

在这里插入图片描述

挺好玩的结论：

质量越大，有量纲的投射距离越大，但是无量纲的投射距离越小。
质量很大之后，有量纲的投射距离会趋于一个定值，但是无量纲的投射距离会趋于0。

因为参考长度为 $m / k$ ，所以质量越大，参考长度越大，这就是上面现象的原因。

至于其他参数的变化，读者可以自行尝试。

完整代码

classdef ProjectileApp < matlab.apps.AppBase
    properties
        % UI Components
        Figure
        MainGrid      % 主网格布局
        LeftGrid      % 左侧网格布局
        RightGrid     % 右侧网格布局
        NonDimAxes       % 无量纲坐标轴
        DimensionalAxes % 有量纲坐标轴
        
        % Tab Group
        TabGroup
        TrajectoryTab
        DimensionalTab
        SurfaceTab
        
        % Input Parameters
        MassEdit
        DragCoeffEdit
        GravityEdit
        VelocityEdit
        AngleEdit
        
        % Simulation Parameters
        Lambda          % 无量纲量 mg/(kV^2)
        LengthScale    % 无量纲长度 m/k
        TimeScale      % 无量纲时间 m/(kV)
        Theta
        
        % Display Labels
        LengthScaleLabel
        TimeScaleLabel
        
        % Plot Data
        TimeData
        XData
        YData
        
        % Value Changed Listeners
        MassListener
        DragCoeffListener
        GravityListener
        VelocityListener
        AngleListener
        
        % 最高点显示
        MaxHeightLabel
        MaxHeightNondimLabel
        
        % 3D Plot
        Surface3DAxes    % 三维图坐标轴
        LambdaRange      % lambda 的范围数组
        ThetaRange      % theta 的范围数组
        RangeData      % 存储射程数据的矩阵
        
        % 轨迹数据存储
        TrajectoryData    % 存储所有轨迹数据的结构体数组
        TrajectoryLines   % 存储所有轨迹线的句柄数组
        LegendEntries     % 存储图例条目
        
        % 删除按钮
        ClearButton
    end
    
    methods
        function app = ProjectileApp
            % 初始化基类
            app = app@matlab.apps.AppBase();
            
            % 创建主窗口
            app.Figure = uifigure('Name', 'Projectile Motion Analysis', ...
                'Position', [100 100 1200 600]);
            movegui(app.Figure, 'center');
            
            % 创建主网格布局 - 2列，左窄右宽
            app.MainGrid = uigridlayout(app.Figure, [1 2]);
            app.MainGrid.ColumnWidth = {'3x', '7x'};
            
            % 创建左侧控制面板网格 - 增加行数以容纳无量纲量显示
            app.LeftGrid = uigridlayout(app.MainGrid, [11 2]);
            app.LeftGrid.RowHeight = {'fit', 'fit', 'fit','fit', 'fit', 'fit', 'fit', 'fit', 'fit', 'fit', 'fit'};
            app.LeftGrid.Padding = [10 10 10 10];
            app.LeftGrid.RowSpacing = 10;
            
            % 创建右侧图表网格
            app.RightGrid = uigridlayout(app.MainGrid, [1 1]);
            app.RightGrid.Padding = [10 10 10 10];
            
            % 创建 Tab Group
            app.TabGroup = uitabgroup(app.RightGrid);
            
            % 创建轨迹 Tab
            app.TrajectoryTab = uitab(app.TabGroup);
            app.TrajectoryTab.Title = 'Non-dimensional';
            
            % 创建轨迹 Tab 的网格布局
            trajectoryGrid = uigridlayout(app.TrajectoryTab, [1 1]);
            trajectoryGrid.Padding = [10 10 10 10];
            
            % 创建有量纲轨迹 Tab
            app.DimensionalTab = uitab(app.TabGroup);
            app.DimensionalTab.Title = 'Dimensional';
            
            % 创建有量纲轨迹 Tab 的网格布局
            dimensionalGrid = uigridlayout(app.DimensionalTab, [1 1]);
            dimensionalGrid.Padding = [10 10 10 10];
            
            % 创建三维图 Tab
            app.SurfaceTab = uitab(app.TabGroup);
            app.SurfaceTab.Title = 'Range Surface';
            
            % 创建三维图 Tab 的网格布局
            surfaceGrid = uigridlayout(app.SurfaceTab, [1 1]);
            surfaceGrid.Padding = [10 10 10 10];
            
            % 创建输入控件
            createInputFields(app);
            
            % 创建轨迹图
            app.NonDimAxes = uiaxes(trajectoryGrid);
            app.NonDimAxes.Layout.Row = 1;
            app.NonDimAxes.Layout.Column = 1;
            title(app.NonDimAxes, 'Non-dimensional Trajectories');
            xlabel(app.NonDimAxes, 'X (Non-dimensional)');
            ylabel(app.NonDimAxes, 'Y (Non-dimensional)');
            grid(app.NonDimAxes, 'on');
            
            % 创建有量纲轨迹图
            app.DimensionalAxes = uiaxes(dimensionalGrid);
            app.DimensionalAxes.Layout.Row = 1;
            app.DimensionalAxes.Layout.Column = 1;
            title(app.DimensionalAxes, 'Dimensional Trajectories');
            xlabel(app.DimensionalAxes, 'x (m)');
            ylabel(app.DimensionalAxes, 'y (m)');
            grid(app.DimensionalAxes, 'on');
            
            % 创建三维图
            app.Surface3DAxes = uiaxes(surfaceGrid);
            app.Surface3DAxes.Layout.Row = 1;
            app.Surface3DAxes.Layout.Column = 1;
            title(app.Surface3DAxes, 'Range Surface');
            xlabel(app.Surface3DAxes, 'λ');
            ylabel(app.Surface3DAxes, 'θ (deg)');
            zlabel(app.Surface3DAxes, 'Range (m)');
            grid(app.Surface3DAxes, 'on');
            view(app.Surface3DAxes, 45, 30);
            
            % 初始化三维图
            calculateRangeSurface(app);
            
            % 初始化无量纲参数
            updateDimensionlessParameters(app);
        end
        
        function createInputFields(app)
            % 质量输入
            lbl = uilabel(app.LeftGrid);
            lbl.Text = 'Mass (kg):';
            lbl.Layout.Row = 1;
            lbl.Layout.Column = 1;
            
            app.MassEdit = uislider(app.LeftGrid);
            app.MassEdit.Limits = [0.1 10];
            app.MassEdit.Value = 1;
            app.MassEdit.Layout.Row = 1;
            app.MassEdit.Layout.Column = 2;
            
            % 阻力系数输入
            lbl = uilabel(app.LeftGrid);
            lbl.Text = 'Drag Coefficient (N·s²/m²):';
            lbl.Layout.Row = 2;
            lbl.Layout.Column = 1;
            
            app.DragCoeffEdit = uislider(app.LeftGrid);
            app.DragCoeffEdit.Limits = [0.01 1];
            app.DragCoeffEdit.Value = 0.1;
            app.DragCoeffEdit.Layout.Row = 2;
            app.DragCoeffEdit.Layout.Column = 2;
            
            % 重力加速度输入
            lbl = uilabel(app.LeftGrid);
            lbl.Text = 'Gravity (m/s²):';
            lbl.Layout.Row = 3;
            lbl.Layout.Column = 1;
            
            app.GravityEdit = uislider(app.LeftGrid);
            app.GravityEdit.Limits = [1 20];
            app.GravityEdit.Value = 9.81;
            app.GravityEdit.Layout.Row = 3;
            app.GravityEdit.Layout.Column = 2;
            
            % 初速度输入
            lbl = uilabel(app.LeftGrid);
            lbl.Text = 'Initial Velocity (m/s):';
            lbl.Layout.Row = 4;
            lbl.Layout.Column = 1;
            
            app.VelocityEdit = uislider(app.LeftGrid);
            app.VelocityEdit.Limits = [1 50];
            app.VelocityEdit.Value = 10;
            app.VelocityEdit.Layout.Row = 4;
            app.VelocityEdit.Layout.Column = 2;
            
            % 角度输入
            lbl = uilabel(app.LeftGrid);
            lbl.Text = 'Angle (degrees):';
            lbl.Layout.Row = 5;
            lbl.Layout.Column = 1;
            
            app.AngleEdit = uislider(app.LeftGrid);
            app.AngleEdit.Limits = [0 90];
            app.AngleEdit.Value = 45;
            app.AngleEdit.Layout.Row = 5;
            app.AngleEdit.Layout.Column = 2;
            
            % 无量纲量显示
            lbl = uilabel(app.LeftGrid);
            lbl.Text = 'Length Scale (m/k):';
            lbl.Layout.Row = 6;
            lbl.Layout.Column = 1;
            
            app.LengthScaleLabel = uilabel(app.LeftGrid);
            app.LengthScaleLabel.Text = '0';
            app.LengthScaleLabel.Layout.Row = 6;
            app.LengthScaleLabel.Layout.Column = 2;
            
            lbl = uilabel(app.LeftGrid);
            lbl.Text = 'Time Scale (m/kV):';
            lbl.Layout.Row = 7;
            lbl.Layout.Column = 1;
            
            app.TimeScaleLabel = uilabel(app.LeftGrid);
            app.TimeScaleLabel.Text = '0';
            app.TimeScaleLabel.Layout.Row = 7;
            app.TimeScaleLabel.Layout.Column = 2;
            
            % 添加最高点显示
            lbl = uilabel(app.LeftGrid);
            lbl.Text = 'Max Height (m):';
            lbl.Layout.Row = 8;
            lbl.Layout.Column = 1;
            
            app.MaxHeightLabel = uilabel(app.LeftGrid);
            app.MaxHeightLabel.Text = '0';
            app.MaxHeightLabel.Layout.Row = 8;
            app.MaxHeightLabel.Layout.Column = 2;
            
            lbl = uilabel(app.LeftGrid);
            lbl.Text = 'Max Height (Non-dim):';
            lbl.Layout.Row = 9;
            lbl.Layout.Column = 1;
            
            app.MaxHeightNondimLabel = uilabel(app.LeftGrid);
            app.MaxHeightNondimLabel.Text = '0';
            app.MaxHeightNondimLabel.Layout.Row = 9;
            app.MaxHeightNondimLabel.Layout.Column = 2;
            
            % 创建按钮网格布局
            buttonGrid = uigridlayout(app.LeftGrid, [2 1]);
            buttonGrid.Layout.Row = 10;
            buttonGrid.Layout.Column = [1 2];
            buttonGrid.RowHeight = {'1x', '1x'};
            buttonGrid.RowSpacing = 5;
            
            % 模拟按钮
            btn = uibutton(buttonGrid);
            btn.Text = 'Simulate';
            btn.ButtonPushedFcn = @app.simulateButtonPushed;
            btn.Layout.Row = 1;
            btn.Layout.Column = 1;
            
            % 清空按钮
            app.ClearButton = uibutton(buttonGrid);
            app.ClearButton.Text = 'Clear All';
            app.ClearButton.ButtonPushedFcn = @app.clearButtonPushed;
            app.ClearButton.Layout.Row = 2;
            app.ClearButton.Layout.Column = 1;
            
            % 添加值改变事件监听
            app.MassListener = listener(app.MassEdit, 'ValueChanged', @(~,~) updateDimensionlessParameters(app));
            app.DragCoeffListener = listener(app.DragCoeffEdit, 'ValueChanged', @(~,~) updateDimensionlessParameters(app));
            app.GravityListener = listener(app.GravityEdit, 'ValueChanged', @(~,~) updateDimensionlessParameters(app));
            app.VelocityListener = listener(app.VelocityEdit, 'ValueChanged', @(~,~) updateDimensionlessParameters(app));
            app.AngleListener = listener(app.AngleEdit, 'ValueChanged', @(~,~) updateDimensionlessParameters(app));
        end
        
        function updateDimensionlessParameters(app)
            % 获取当前参数
            m = app.MassEdit.Value;
            k = app.DragCoeffEdit.Value;
            g = app.GravityEdit.Value;
            V = app.VelocityEdit.Value;
            theta_deg = app.AngleEdit.Value;
            
            % 计算无量纲参数
            app.LengthScale = m/k;
            app.TimeScale = m/(k*V);
            app.Lambda = m*g/(k*V^2);
            app.Theta = theta_deg * pi/180;
            
            % 更新显示
            app.LengthScaleLabel.Text = sprintf('%.2f m', app.LengthScale);
            app.TimeScaleLabel.Text = sprintf('%.2f s', app.TimeScale);
            
            % 检查是否已存在相同无量纲参数的轨迹
            if ~isempty(app.TrajectoryData)
                existingParams = [app.TrajectoryData.params];
                for i = 1:length(existingParams)
                    % 计算已存在轨迹的无量纲参数
                    existingLambda = existingParams(i).m * existingParams(i).g / ...
                        (existingParams(i).k * existingParams(i).V^2);
                    existingTheta = existingParams(i).theta * pi/180;
                    
                    % 使用无量纲参数进行比较
                    if abs(existingLambda - app.Lambda) < 1e-6 && ...
                            abs(existingTheta - app.Theta) < 1e-6
                        % 如果找到相同无量纲参数，直接返回
                        return;
                    end
                end
            end
            
            % 求解ODE
            [app.TimeData, app.XData, app.YData] = solveDimensionlessODE(app);
            
            % 计算最高点
            [maxHeight, maxIdx] = max(app.YData);
            maxHeightDim = maxHeight * app.LengthScale;
            
            % 更新最高点显示
            app.MaxHeightLabel.Text = sprintf('%.2f m', maxHeightDim);
            app.MaxHeightNondimLabel.Text = sprintf('%.2f', maxHeight);
            
            % 创建新的轨迹线 - 使用无量纲坐标
            hold(app.NonDimAxes, 'on');
            newLine = plot(app.NonDimAxes, app.XData, app.YData, 'LineWidth', 2);
            hold(app.NonDimAxes, 'off');
            
            % 创建新的轨迹线 - 使用有量纲坐标
            hold(app.DimensionalAxes, 'on');
            newDimensionalLine = plot(app.DimensionalAxes, app.XData * app.LengthScale, app.YData * app.LengthScale, 'LineWidth', 2);
            hold(app.DimensionalAxes, 'off');
            
            % 存储轨迹数据 - 同时存储有量纲和无量纲坐标
            newData = struct('x', app.XData * app.LengthScale, 'y', app.YData * app.LengthScale, ...
                'X', app.XData, 'Y', app.YData, ...
                'params', struct('m', m, 'k', k, 'g', g, 'V', V, 'theta', theta_deg));
            app.TrajectoryData = [app.TrajectoryData, newData];
            app.TrajectoryLines = [app.TrajectoryLines, newLine, newDimensionalLine];
            
            % 更新图例
            legendText = sprintf('λ=%.2f, θ=%.1f°', app.Lambda, theta_deg);
            app.LegendEntries = [app.LegendEntries, {legendText}];
            
            % 在最高点添加文本标注
            [maxY, maxIdx] = max(app.YData);
            text(app.NonDimAxes, app.XData(maxIdx), maxY, legendText, ...
                'VerticalAlignment', 'bottom', ...
                'HorizontalAlignment', 'center');
            
            % 在有量纲图中添加最高点标注
            text(app.DimensionalAxes, app.XData(maxIdx) * app.LengthScale, maxY * app.LengthScale, legendText, ...
                'VerticalAlignment', 'bottom', ...
                'HorizontalAlignment', 'center');
            
            % 将图例移到坐标轴外部
            legend(app.NonDimAxes, app.LegendEntries, 'Location', 'eastoutside');
            legend(app.DimensionalAxes, app.LegendEntries, 'Location', 'eastoutside');
            
            % 设置坐标轴范围为自动
            axis(app.NonDimAxes, 'auto');
            axis(app.DimensionalAxes, 'auto');
            
            % 更新标题
            title(app.NonDimAxes, sprintf('Non-dimensional Trajectories (λ=%.2f, θ=%.1f°)', app.Lambda, theta_deg));
            title(app.DimensionalAxes, sprintf('Dimensional Trajectories (λ=%.2f, θ=%.1f°)', app.Lambda, theta_deg));
            
            % 更新三维图
            calculateRangeSurface(app);
        end
        
        function simulateButtonPushed(app, ~, ~)
            % 直接调用updateDimensionlessParameters来更新轨迹
            updateDimensionlessParameters(app);
        end
        
        function [T, X, Y] = solveDimensionlessODE(app)
            % 初始条件
            X0 = 0;
            Y0 = 0;
            Xp0 = cos(app.Theta);
            Yp0 = sin(app.Theta);
            
            % 求解ODE直到Y回到初始高度
            y0 = [X0; Y0; Xp0; Yp0];
            options = odeset('Events', @(t,y) heightEvent(t,y,Y0), ...
                'RelTol', 1e-8, ...  % 相对误差
                'AbsTol', 1e-8, ...  % 绝对误差
                'MaxStep', 0.01);    % 最大步长
            [T, Y] = ode45(@(t,y) dimensionlessODE(t, y, app.Lambda), [0 100], y0, options);
            
            X = Y(:,1);
            Y = Y(:,2);
        end
        
        function calculateRangeSurface(app)
            % 创建lambda和theta的网格
            app.LambdaRange = linspace(0.01, 100, 20);
            app.ThetaRange = linspace(1, 90, 20);  % 从1度开始，而不是0度
            [Lambda, Theta] = meshgrid(app.LambdaRange, app.ThetaRange);
            app.RangeData = zeros(size(Lambda));
            
            % 计算每个点的最终x距离
            for i = 1:length(app.ThetaRange)
                for j = 1:length(app.LambdaRange)
                    lambda = Lambda(i,j);
                    theta = Theta(i,j) * pi/180;
                    
                    % 求解ODE直到Y回到初始高度
                    y0 = [0; 0; cos(theta); sin(theta)];
                    options = odeset('Events', @(t,y) heightEvent(t,y,0), ...
                        'RelTol', 1e-8, ...
                        'AbsTol', 1e-8, ...
                        'MaxStep', 0.01);
                    [~, Y] = ode45(@(t,y) dimensionlessODE(t, y, lambda), [0 100], y0, options);
                    
                    % 检查是否有解
                    if ~isempty(Y)
                        
                        % 获取最终x距离（转换为有量纲）
                        % Y(end,1)是物体落地时的x坐标，Y(end,2)是y坐标（应该接近0）
                        if abs(Y(end,2)) < 1e-6  % 确保物体确实落地
                            app.RangeData(i,j) = Y(end,1) * lambda;
                        else
                            % 如果物体没有正确落地，设置为NaN
                            app.RangeData(i,j) = NaN;
                        end
                    else
                        % 如果没有解，设置为NaN
                        app.RangeData(i,j) = NaN;
                    end
                end
            end
            
            % 绘制三维图
            surf(app.Surface3DAxes, Lambda, Theta, app.RangeData);
            colorbar(app.Surface3DAxes);
            shading(app.Surface3DAxes, 'interp');
            
            % 更新标签
            title(app.Surface3DAxes, 'Range Surface');
            xlabel(app.Surface3DAxes, 'λ');
            ylabel(app.Surface3DAxes, 'θ (deg)');
            zlabel(app.Surface3DAxes, 'Range (m)');
            
            % 在当前参数位置添加标记
            % hold(app.Surface3DAxes, 'on');
            % plot3(app.Surface3DAxes, app.Lambda, app.Theta*180/pi, ...
            %     app.XData(end) * app.LengthScale, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
            % hold(app.Surface3DAxes, 'off');
        end
        
        function clearButtonPushed(app, ~, ~)
            % 清除所有轨迹线
            if ~isempty(app.TrajectoryLines)
                delete(app.TrajectoryLines);
                app.TrajectoryLines = [];
            end
            
            % 清除所有轨迹数据
            app.TrajectoryData = [];
            
            % 清除图例
            app.LegendEntries = {};
            legend(app.NonDimAxes, 'off');
            legend(app.DimensionalAxes, 'off');
            
            % 清除坐标轴上的文本标注
            for ax = [app.NonDimAxes, app.DimensionalAxes]
                children = get(ax, 'Children');
                for i = 1:length(children)
                    if isa(children(i), 'matlab.graphics.primitive.Text')
                        delete(children(i));
                    end
                end
            end
            
            % 重置坐标轴
            axis(app.NonDimAxes, 'auto');
            axis(app.DimensionalAxes, 'auto');
            
            % 更新标题
            title(app.NonDimAxes, 'Non-dimensional Trajectories');
            xlabel(app.NonDimAxes, 'X (Non-dimensional)');
            ylabel(app.NonDimAxes, 'Y (Non-dimensional)');
            
            title(app.DimensionalAxes, 'Dimensional Trajectories');
            xlabel(app.DimensionalAxes, 'x (m)');
            ylabel(app.DimensionalAxes, 'y (m)');
            
            % 更新三维图
            calculateRangeSurface(app);
        end
    end
end

function [value, isterminal, direction] = heightEvent(t, y, Y0)
% 当Y回到初始高度时停止
value = y(2) - Y0;      % 当前高度与初始高度的差
isterminal = 1;         % 停止积分
direction = -1;         % 只在下降穿过初始高度时停止
end

function dydt = dimensionlessODE(t, y, lambda)
% y = [X; Y; X'; Y']
X = y(1);
Y = y(2);
Xp = y(3);
Yp = y(4);

% Calculate velocity magnitude
V = sqrt(Xp^2 + Yp^2);

% ODE system
dydt = zeros(4,1);
dydt(1) = Xp;
dydt(2) = Yp;
dydt(3) = -V * Xp;
dydt(4) = -V * Yp - lambda;
end