Dimensional Analysis量纲分析入门

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

幂律（Power Law）

在物理学中，幂律（Power Law）是指量之间的幂次关系，通常表示为：
$$y=k{x}_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_{1}f(\alpha x_2)f^{\prime }(\alpha x_1)=x_{2}f(\alpha x_1)f^{\prime }(\alpha x_2)$$

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

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

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

$$\int \frac{f^{\prime }(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$）, 单位：千克（$kg$）
• 时间（$T$）, 单位：秒（$s$）
• 电流（$I$）, 单位：安培（$A$）
• 热力学温度（$\Theta$）, 单位：开尔文（$K$）
• 物质的量（$N$）, 单位：摩尔（$mol$）
• 发光强度（$J$）, 单位：坎德拉（$cd$）

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

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

$$y=C{a}^{\alpha }{b}^{\beta }{c}^{\gamma }\cdots$$

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

$$[y]=[a{]}^{\alpha }[b{]}^{\beta }[c{]}^{\gamma }\cdots$$

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

$$\begin{array}{rl}& [y]=L^{{\alpha }_y}{M}^{{\beta }_y}{T}^{{\gamma }_y}{I}^{{\delta }_y}{\Theta }^{{ϵ}_y}{N}^{{\zeta }_y}{J}^{{\eta }_y}\\ & [a]=L^{{\alpha }_a}{M}^{{\beta }_a}{T}^{{\gamma }_a}{I}^{{\delta }_a}{\Theta }^{{ϵ}_a}{N}^{{\zeta }_a}{J}^{{\eta }_a}\\ & [b]=L^{{\alpha }_b}{M}^{{\beta }_b}{T}^{{\gamma }_b}{I}^{{\delta }_b}{\Theta }^{{ϵ}_b}{N}^{{\zeta }_b}{J}^{{\eta }_b}\\ & [c]=L^{{\alpha }_c}{M}^{{\beta }_c}{T}^{{\gamma }_c}{I}^{{\delta }_c}{\Theta }^{{ϵ}_c}{N}^{{\zeta }_c}{J}^{{\eta }_c}\\ & \cdots \end{array}$$

这样可以得到一组方程：

$$\begin{array}{rl}& {\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 \\ & {ϵ}_y={ϵ}_a\cdot \alpha +{ϵ}_b\cdot \beta +{ϵ}_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{array}$$

举个例子

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

我们写成：

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

量纲分析有：

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

因此有：

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

$$\{\begin{array}{rl}& \alpha +\gamma =0\\ & \beta =0\\ & -2\gamma =-1\end{array}$$

解得：

$$\{\begin{array}{rl}& \alpha =-\frac{1}{2}\\ & \beta =0\\ & \gamma =\frac{1}{2}\end{array}$$

所以有：

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

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

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

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

周期为

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

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

进一步的讨论

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

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

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

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

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

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

$$\{\begin{array}{rl}& \alpha +\gamma =0\\ & \beta -\gamma +\delta =1\\ & -2\beta -\delta =0\end{array}$$

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

$$\{\begin{array}{rl}& \alpha =1\\ & \beta =0\\ & \gamma =-1\\ & \delta =0\end{array}$$

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

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

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

$$\{\begin{array}{rl}& \alpha +\gamma =0\\ & \beta -\gamma +\delta =0\\ & -2\beta -\delta =0\end{array}$$

解得：

$$\{\begin{array}{rl}& \beta =\alpha \\ & \gamma =-\alpha \\ & \delta =-2\alpha \end{array}$$

所以有：

$$\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})$$

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

例如，

$$\{\begin{array}{rl}& \alpha +\gamma =0\\ & \beta -\gamma +\delta =1\\ & -2\beta -\delta =0\end{array}$$

的解还可能是：

$$\{\begin{array}{rl}& \alpha =0\\ & \beta =-1\\ & \gamma =0\\ & \delta =2\end{array}$$

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

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

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

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

观察发现，

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

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

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

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

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

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

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

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

或者

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

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

$$h=\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$。

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

展开得到：

$$\begin{array}{rl}& 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{array}$$

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

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

$$X=\frac{x}{m/k},Y=\frac{y}{m/k},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}^{\prime }$$

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

最终可以得到：

$$\begin{array}{rl}& k{V}^2{X}^{\prime \prime }=-k\sqrt{V^2{X}^{\prime 2}+V^2{Y}^{\prime 2}}V{X}^{\prime }\\ & k{V}^2{Y}^{\prime \prime }=-k\sqrt{V^2{X}^{\prime 2}+V^2{Y}^{\prime 2}}V{Y}^{\prime }-mg\end{array}$$

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

$$\begin{array}{rl}& X^{\prime \prime }=-\sqrt{X^{\prime 2}+Y^{\prime 2}}{X}^{\prime }\\ & Y^{\prime \prime }=-\sqrt{X^{\prime 2}+Y^{\prime 2}}{Y}^{\prime }-\lambda \end{array}$$

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

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

程序实现

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

function [value, isterminal, direction] = heightEvent(t, y, Y0)
% 当Y回到初始高度时停止
value = y(2) - Y0;      % 当前高度与初始高度的差
isterminal = 1;         % 停止积分
direction = -1;         % 只在下降穿过初始高度时停止
endfunction 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.AppBaseproperties% UI ComponentsFigureMainGrid      % 主网格布局LeftGrid      % 左侧网格布局RightGrid     % 右侧网格布局NonDimAxes       % 无量纲坐标轴DimensionalAxes % 有量纲坐标轴% Tab GroupTabGroupTrajectoryTabDimensionalTabSurfaceTab% Input ParametersMassEditDragCoeffEditGravityEditVelocityEditAngleEdit% Simulation ParametersLambda          % 无量纲量 mg/(kV^2)LengthScale    % 无量纲长度 m/kTimeScale      % 无量纲时间 m/(kV)Theta% Display LabelsLengthScaleLabelTimeScaleLabel% Plot DataTimeDataXDataYData% Value Changed ListenersMassListenerDragCoeffListenerGravityListenerVelocityListenerAngleListener% 最高点显示MaxHeightLabelMaxHeightNondimLabel% 3D PlotSurface3DAxes    % 三维图坐标轴LambdaRange      % lambda 的范围数组ThetaRange      % theta 的范围数组RangeData      % 存储射程数据的矩阵% 轨迹数据存储TrajectoryData    % 存储所有轨迹数据的结构体数组TrajectoryLines   % 存储所有轨迹线的句柄数组LegendEntries     % 存储图例条目% 删除按钮ClearButtonendmethodsfunction 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 Groupapp.TabGroup = uitabgroup(app.RightGrid);% 创建轨迹 Tabapp.TrajectoryTab = uitab(app.TabGroup);app.TrajectoryTab.Title = 'Non-dimensional';% 创建轨迹 Tab 的网格布局trajectoryGrid = uigridlayout(app.TrajectoryTab, [1 1]);trajectoryGrid.Padding = [10 10 10 10];% 创建有量纲轨迹 Tabapp.DimensionalTab = uitab(app.TabGroup);app.DimensionalTab.Title = 'Dimensional';% 创建有量纲轨迹 Tab 的网格布局dimensionalGrid = uigridlayout(app.DimensionalTab, [1 1]);dimensionalGrid.Padding = [10 10 10 10];% 创建三维图 Tabapp.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);endfunction 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));endfunction 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;endendend% 求解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);endfunction simulateButtonPushed(app, ~, ~)% 直接调用updateDimensionlessParameters来更新轨迹updateDimensionlessParameters(app);endfunction [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);endfunction 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% 如果物体没有正确落地，设置为NaNapp.RangeData(i,j) = NaN;endelse% 如果没有解，设置为NaNapp.RangeData(i,j) = NaN;endendend% 绘制三维图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');endfunction 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));endendend% 重置坐标轴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);endend
endfunction [value, isterminal, direction] = heightEvent(t, y, Y0)
% 当Y回到初始高度时停止
value = y(2) - Y0;      % 当前高度与初始高度的差
isterminal = 1;         % 停止积分
direction = -1;         % 只在下降穿过初始高度时停止
endfunction 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

总结

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