专栏文章

使用 RK 法求解变分问题所对应的微分方程：P13210 [GCJ 2016 Finals] Radioactive Islands

kk43a1_sad
P13210题解参与者 3已保存评论 2

文章操作

快速查看文章及其快照的属性，并进行相关操作。

当前评论
2 条
当前快照
1 份
快照标识符
@mio77o1a
此快照首次捕获于
2025/12/02 14:29
3 个月前
此快照最后确认于
2025/12/02 14:29
3 个月前
查看原文

洛谷(P13210) Radioactive Islands ：使用 RK 法求解变分问题所对应的微分方程

相比于其他方法，此方法计算速度可以较快，使用最后给出的代码，对于此题目的数据集，总耗时约 37ms37\mathrm{ms}。
起点终点已确定，题目要求中，要最小化的泛函是： J=∫AB(1+∑i=0N−1Di−2)ds,Di=x2+(y−Ci)2J=\int_A^{B}\left(1+\sum_{i=0}^{N-1} D_i^{-2}\right)\mathrm{d}s,\quad D_i=\sqrt{x^2+(y-C_i)^2} 我们假设 CiC_i 已从小到大排序，是辐射源的位置，即 i<j⇒Ci<Cji<j\Rightarrow C_i<C_j 其最小值必然满足变分为 00 即 δJ=0\delta J=0，我们可以用变分为 00 这个条件排除许多路径。类比费马原理，即光程变分为 00： δS=δ∫ABnds=0\delta S=\delta \int_{A}^B n\mathrm{d}s=0 可以看到，如果我们令 n(x,y)=1+∑i=0N−1Di−2n(x,y)=1+\sum_{i=0}^{N-1} D_i^{-2} 那么费马原理就和原问题所需最小化函数的变分为 00 形式上相同。于是可以知道，原问题中所受辐射最小路径必然和某一条在折射率为 n(x,y)n(x,y) 的空间中通过 AB 点的光线重合。于是我们可以尝试以 A 为起点，在初始方向与 xx 正方向夹角在 π2\frac{\pi}{2} 范围内的所有光线进行搜索，看看有没有通过 B 点的。可以预见的，这样的光线应该只有有限根通过 B 点，那么进而直接在这有限的几条光线中求出所受辐射最小路径即可。

1. 光线方程推导

首先我们规定符号： y′=dydx,y′′=d2ydx2y'=\frac{\mathrm{d}y}{\mathrm{d}x},\quad y''=\frac{\mathrm{d}^2y}{\mathrm{d}x^2}
从一般的光线方程出发（此方程也可以从费马原理开始使用欧拉-拉格朗日方程得到，要注意的是有半完整约束 drds⋅drds=(dxds)2+(dyds)2=1\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}s}\cdot \frac{\mathrm{d}\mathbf{r}}{\mathrm{d}s}=\left(\frac{\mathrm{d}x}{\mathrm{d}s}\right)^2+\left(\frac{\mathrm{d}y}{\mathrm{d}s}\right)^2=1 ，此半完整约束的处理办法是使用带拉格朗日乘子的欧拉-拉格朗日方程，通过一些化简与边界条件可以得到拉格朗日乘子 2λ=n2\lambda=n ）： dds(ndrds)=∇n\frac{\mathrm{d}}{\mathrm{d}s}\left(n\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}s}\right)=\nabla n 其中 ss 为光线的曲线长度，r=(x(s),y(s))\mathbf{r}=(x(s),y(s)) 是光线的参数方程，以 ss 为参数，可以知道 drds\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}s} 即为光线的单位切向量。进行一些计算可以得到以 xx 为自变量的光线微分方程： y′′=1+y′2n(∂n∂y−y′∂n∂x)y''=\frac{1+y'^2}{n}\left(\frac{\partial n}{\partial y}-y'\frac{\partial n}{\partial x}\right)
注：实际上以弧长 ss 为自变量可能是更为自然的选择，对应的方程在最后的（7.(2) 从光线方程到以 xx 为自变量的微分方程），可能会让 RK 法步数更少，以及奇点附近表现更好，但是我一开始使用的是 xx，所以就将就了，以弧长 ss 为自变量留给读者尝试。
折射率的偏导也不难求： ∂n∂x=−∑i=0N−12x(x2+(y−Ci)2)2\frac{\partial n}{\partial x}=-\sum_{i=0}^{N-1} \frac{2x}{(x^2+(y-C_i)^2)^2} ∂n∂y=−∑i=0N−12(y−Ci)(x2+(y−Ci)2)2\frac{\partial n}{\partial y}=-\sum_{i=0}^{N-1} \frac{2(y-C_i)}{(x^2+(y-C_i)^2)^2} 高阶方程皆可降次，只需把 y′y' 看作未知函数，且遵守微分方程： ddxy′=y′′\frac{\mathrm{d}}{\mathrm{d}x}y'=y'' ddxy=y′\frac{\mathrm{d}}{\mathrm{d}x}y=y' 我们在计算路径的过程中，可以顺便算出所受辐射量，只需要把 JJ 看成未知函数，添加一个微分方程： dJdx=n1+y′2\frac{\mathrm{d}J}{\mathrm{d}x}=n\sqrt{1+y'^2} 最终的微分方程组为： {dydx=y′,y(−10)=Ady′dx=1+y′2n(∂n∂y−y′∂n∂x),y′(−10)=tan⁡θdJdx=n1+y′2,J(−10)=0\begin{cases} \frac{\mathrm{d}y}{\mathrm{d}x}=y' ,& y(-10)=A \\ \frac{\mathrm{d}y'}{\mathrm{d}x}=\frac{1+y'^2}{n}\left(\frac{\partial n}{\partial y}-y'\frac{\partial n}{\partial x}\right),& y'(-10)=\tan{\theta} \\ \frac{\mathrm{d}J}{\mathrm{d}x}=n\sqrt{1+y'^2},& J(-10)=0 \end{cases}

2. 显式 Runge-Kutta 法解微分方程

我这里用的是 Bogacki 和 Shampine 的 5(4)5(4) 格式，具体的系数表可以在 mathematica 中输入：
CPP
NDSolve`EmbeddedExplicitRungeKuttaCoefficients[5, Infinity]
这是一种具有嵌入对的 RK 法，也就是计算一步可以得到一个 44 阶解和一个 55 阶解，55 阶解用于得到下一步的值，55 阶减去 44 阶可以得到一个对误差的估计，使用这个误差的估计，并利用 pp 阶方法的单步误差正比于步长 hh 的 hp+1h^{p+1} 的特点，从设定的误差限得到下一步的步长： hn+1=hn(Tol∣∣err∣∣)1/5h_{n+1}=h_n\left(\frac{\mathrm{Tol}}{||\mathrm{err}||}\right)^{1/5} 以此在保证误差的情况下尽量减少迭代步数。

3. 微分方程的停止条件

首先是当 x=10x=10 时停止，这是当然的。
另外我们需要计算 x=0x=0 时 y(0)y(0) 的值，这帮助我们分类区间，具体将在之后讲解。
另外有几种情况也是需要停止的，首先我们对最短路径的所受辐射上界作一估计：
Jmax=(1+0.01N)(min⁡(20−A−B,20+A+B)+10π)J_{max}=(1 + 0.01N )(\min(20 - A - B, 20 + A + B) + 10\pi ) 这相当于沿着 x=10,x=−10x=10,x=-10 两条线走到 y=±10y=\pm 10，绕一个半圆走过去，路上的折射率一定小于 (1+0.01N)(1 + 0.01N )。

(1) yy 超限

那么当 yy 过于大，以至于 ∣y∣−10>Jmax/2|y|-10>J_{max}/2 时显然不是所受辐射最小路径，停止计算。

(2) JJ 超限

另外由于我们在计算路径的过程中可以计算已受辐射量，那么可以对现在路径受到辐射的下界进行一个估计： J+(x−10)2+(y−B)2<JmaxJ+\sqrt{(x-10)^2+(y-B)^2}<J_{max} 如果上式不成立且 x>0x>0，那么可以停止计算，因为然不是所受辐射最小路径。
另外两种停止条件没那么直观：

(3) y′y' 发散

第一种情况，在初始角度随意的情况下，光线可能往回拐，这会导致 y′y' 在有限的x内发散到无穷。但是就如题目的 pdf 附件中说的，我们往回拐的一定不是所受辐射最小路径，因为可以沿着 x=cosntx=cosnt 连接拐回去那一点，这样得到的路径更短。
判断这种情况的办法是计算当前的光线与 yy 轴的夹角 θ\theta，再计算当前的曲率半径： R=∣(1+y′2)32y′′∣R=|\frac{(1+y'^2)^{\frac{3}{2}}}{y''}|，如果 RθR\theta 足够小，意味着光线走过路径很短，此时可近似为一圆形，那么在 hh 的 xx 方向步长内 y′y' 发散的条件即是 R(1−cos⁡θ)<hR(1-\cos{\theta})<h 即圆的最右侧在距离当前位置不到 hh 处，此时也应该停止。

(4) 落入辐射源

第二种情况，在初始角度随意的情况下，根据微分方程演化的光线有可能撞上某个辐射源所在的点，这时所受辐射量会发散到无穷，这显然不是所受辐射最小路径。
我们从极坐标的角度来了解这件事，当我们距离某个辐射源足够近的时候，Di−2D_i^{-2} 将迅速发散到无穷，这时由于距离足够近，所以其他辐射源的导致的辐射变化可忽略，可假设 n≈H+Di−2n\approx H+D_i^{-2}，其中 HH 为某个常数，我们换到以 CiC_i 为原点的极坐标系中，可知 n(r)=H+r−2n(r)=H+r^{-2}，折射率只与 rr 有关，是旋转对称的，在这种情况下，有 rn(r)sin⁡θ=constrn(r)\sin{\theta}=const 其中 sin⁡θ=rdϕr2dϕ2+dr2=11+(1rdrdϕ)2\sin{\theta}=\frac{r\mathrm{d}\phi}{\sqrt{r^2\mathrm{d}\phi^2+\mathrm{d}r^2}}=\frac{1}{\sqrt{1+(\frac{1}{r}\frac{\mathrm{d}r}{\mathrm{d}\phi})^2}} 是光线与极坐标系径向的夹角。这个式子我们不做证明，但有两种理解方式，一种是直接对光线进行微元分析可证，另一种理解方式是这个式子代表了光子的角动量守恒，因为系统旋转对称。
进一步的，代入 sin⁡θ\sin{\theta} 与初始条件，上式可化为
drdϕ=−(Hr2+1Hr02+1r0sin⁡θ0)2−r2\frac{\mathrm{d}r}{\mathrm{d}\phi}=-\sqrt{(\frac{Hr^2+1}{Hr_0^2+1}\frac{r_0}{\sin{\theta_0}})^2-r^2} 负号是因为我们假设光线向内运动，在 0<r<r00<r<r_0 时对等号右侧的上限进行估计可知 (drdϕ)2≥(1(H+r0−2)r0sin⁡θ0)2−r02\left(\frac{\mathrm{d}r}{\mathrm{d}\phi}\right)^2\ge (\frac{1}{(H+r_0^{-2})r_0\sin{\theta_0}})^2-r_0^2 如果方程右侧大于 00，则代表 drdϕ\frac{\mathrm{d}r}{\mathrm{d}\phi} 小于某个有限负数，这意味着 rr 将在有限的 ϕ\phi 变化内归 00，即撞上辐射源。于是当足够接近某个辐射源时，方程右侧的式子可作为撞上辐射源的判定依据。

4. 将问题转化为求解非线性方程

直到现在我们讨论的都是在给定初始点 A 与初始光线角度 θ\theta 的情况下，如何求出光线轨迹。而我们最终的目的是找到某几个 θ\theta，使得射出的光线正好经过 B 点，我们可以把 x=10x=10 作为截面，设 y(θ,x)y(\theta,x) 是光线轨迹，已用 RK 法求出，那么当光线没有触发停止条件并到达 x=10x=10 时，函数 g(θ)=y(θ,10)−Bg(\theta)=y(\theta,10)-B 的零点即为我们所求的使得光线正好经过 B 点的 θ\theta，假设我们有到达 x=10x=10 的两点，且 g(θ1)<0,g(θ2)>0g(\theta_1)<0,g(\theta_2)>0，我们就可以用丰富的非线性方程求根办法进行求根，这里我使用的是 Brent 法，在网上扒的。

5. g(θ)g(\theta) 的性质

DBL_MAX 在此代指浮点数能够表示的最大有限值。
我们之前提到了很多终止条件，他们都可能使得光线无法到达 x=10x=10，导致 g(θ)g(\theta) 值无意义，但有时这些终止条件有迹可循，请注意，以下是我通过遍历某些例子的 θ\theta 得到的经验结论，没有严格证明，随着 θ\theta 单调增加，有循环的过程：
⋯\cdots ->落入 Ci−1C_{i-1} 点-> y′y' 发散到负无穷-> g(θ)<0g(\theta)<0 -> g(θ)=0g(\theta)=0 -> g(θ)>0g(\theta)>0 -> y′y' 发散到正无穷->落入 CiC_{i} 点-> ⋯\cdots
且在落入 Ci−1C_{i-1} 点与落入 CiC_{i} 点间，Ci−1≤y(0)≤CiC_{i-1}\le y(0)\le C_{i}。
也就是说，如果规定 y′y' 发散到负无穷时 g(θ)=−g(\theta)=- DBL_MAX，y′y' 发散到正无穷时 g(θ)=g(\theta)= DBL_MAX，落入某一点时 g(θ)=g(\theta)= DBL_MAX，那么无论如何 g(θ),y(θ,0)g(\theta),y(\theta,0) 都有值，且被分割为 N+1N+1 个区间，每个区间上 g(θ)g(\theta) 都单调，在每个区间上有 Ci−1<y(θ,0)≤CiC_{i-1}<y(\theta,0)\le C_{i}，规定 C−1=−∞,CN=∞C_{-1}=-\infty,C_{N}=\infty。
综上所述，我们可以通过 y(θ,0)y(\theta,0) 判断 θ\theta 属于哪个区间，且每个区间内 g(θ)g(\theta) 都单调，可用二分法求解。

6. 二分法预处理

既然有 N+1N+1 个区间就意味着有 N+1N+1 个解，且我们可以知道 θ\theta 属于哪个区间，定义区间编号 ii 满足 Ci−1<y(θi,0)≤CiC_{i-1}<y(\theta_i,0)\le C_{i}，可以把所有区间初始化为 [−π2,π2][-\frac{\pi}{2},\frac{\pi}{2}]，定义这两点 g(−π2)=−g(-\frac{\pi}{2})=- DBL_MAX， g(π2)=g(\frac{\pi}{2})= DBL_MAX，端点所属区间编号为 0,N0,N。

(1) 二分法使区间左右端点的区间编号与区间本身编号重合

对第一个区间进行二分法，得到区间中点 θmid\theta_{mid} 与其所属区间 ii，之后对于所有区间编号小于 ii 的区间右侧进行更新，对所有区间编号大于 ii 的区间左侧进行更新，对第 ii 个区间左右值进行更新，由于每个区间内 g(θ)g(\theta) 都单调增，如果 g(θ)>0g(\theta)>0 更新右值，否则更新左值。这样使得所有的区间不断变小，接着继续对第一个区间进行二分，直到第一个区间的区间编号都为 00。
接着对第二个区间进行相同的处理，但此时不再更新第一个区间。这样处理完所有的区间后，所有的区间内都有一个解，且端点的区间编号都与其编号相同。

(2) 二分法使区间左右端点满足 Brent 法条件

接着对每个区间进行二分法，直到区间两端的 ∣g(θ)∣|g(\theta)| 不等于 DBL_MAX 为止，此时已满足 Brent 法的条件，即可直接求解出第 ii 个区间穿过 B 点的光线对应的 θi\theta_i，在使用 RK 法计算 g(θi)g(\theta_i) 已顺便计算总辐射量 JiJ_i，于是直接求所有 JiJ_i 最小值即可。于是得到所受最小辐射。

7. 繁杂的数学推导

题解中我跳过了一些证明，如果你有兴趣，如下是详细证明。

(1) 从费马原理到光线方程

我们从一般的泛函开始
J[y(x)]=∫x0x1L(x,y(x),y′(x))dx,y(x)=(y1(x),y2(x),…yn(x))J[\mathbf{y}(x)]=\int_{x_0}^{x_1}L(x,\mathbf{y}(x),\mathbf{y}'(x))\mathrm{d}x,\quad \mathbf{y}(x)=(y_1(x),y_2(x),\ldots y_n(x))
当其有 mm 个约束 gi(x,y(x),y′(x))=0g_i(x,\mathbf{y}(x),\mathbf{y}'(x))=0 时，要求变分 δJ=0\delta J=0，令函数 y(x)\mathbf{y}(x) 发生无穷小变化（终点的 x1x_1 不固定）：
y(x),x∈[x0,x1]→y(x)+δy(x),x∈[x0,x1+δx1]\mathbf{y}(x),x\in [x_0,x_1] \to \mathbf{y}(x)+\delta\mathbf{y}(x),x\in [x_0,x_1+\delta x_1]
规定符号：
δy1=δy(x1+δx1)+y(x1+δx1)−y(x1),δy0=δy(x0)\delta\mathbf{y}_1=\delta\mathbf{y}(x_1+\delta x_1)+\mathbf{y}(x_1+\delta x_1)-\mathbf{y}(x_1),\quad \delta\mathbf{y}_0=\delta\mathbf{y}(x_0)
要注意，由于 y′(x)\mathbf{y}'(x) 存在，δx1,δy1,δy(x1)\delta x_1,\delta\mathbf{y}_1,\delta\mathbf{y}(x_1) 不互相独立，它们的线性主部有关系：
δy1=δy(x1+δx1)+y(x1+δx1)−y(x1)=δy(x1)+y′(x1)δx1\delta\mathbf{y}_1=\delta\mathbf{y}(x_1+\delta x_1)+\mathbf{y}(x_1+\delta x_1)-\mathbf{y}(x_1)=\delta\mathbf{y}(x_1)+\mathbf{y}'(x_1)\delta x_1
由于存在约束，则对于变分 δJ=0\delta J=0 时，存在待定函数 λi(x)\lambda_i(x) ，使得我们构造辅助泛函的变分 δJ⋆[y(x)]=0\delta J^{\star}[\mathbf{y}(x)]=0 （如果你把函数看成无穷维向量，那么这一步等同于多元微积分中的拉格朗日乘子法）：
J⋆[y(x)]=∫x0x1(L+∑i=1mλi(x)gi)dx=∫x0x1HdxJ^{\star}[\mathbf{y}(x)]=\int_{x_0}^{x_1}\left(L+\sum_{i=1}^{m}\lambda_i(x)g_i\right)\mathrm{d}x=\int_{x_0}^{x_1}H\mathrm{d}x
规定符号：
Hyi=∂H∂yi,Hyi′=∂H∂yi′,H∣x=t=H(t,y(t)+δy(t),y′(t)+δy′(t))H_{y_i}=\frac{\partial H}{\partial y_i},\quad H_{y'_i}=\frac{\partial H}{\partial y'_i},\quad H|_{x=t}=H(t,\mathbf{y}(t)+\delta\mathbf{y}(t),\mathbf{y}'(t)+\delta\mathbf{y}'(t))
接下来，我们对辅助泛函进行变分，令函数 y(x)\mathbf{y}(x) 发生无穷小变化，辅助泛函的增量为（只保留线性主部）：
δJ⋆=∫x0x1+δx1H(x,y(x)+δy(x),y′(x)+δy′(x))dx−∫x0x1H(x,y(x),y′(x))dx=∫x1x1+δx1H(x,y(x)+δy(x),y′(x)+δy′(x))dx+∫x0x1[H(x,y(x)+δy(x),y′(x)+δy′(x))−H(x,y(x),y′(x))]dx=H∣x=x1δx1+∑i=1n(∫x0x1[Hyiδyi(x)+Hyi′δyi′(x)]dx)=H∣x=x1δx1+∑i=1n(∫x0x1Hyiδyi(x)dx+∫x0x1Hyi′δyi′(x)dx)=H∣x=x1δx1+∑i=1n(∫x0x1Hyiδyi(x)dx+∫x0x1Hyi′dδyi(x)dxdx)=H∣x=x1δx1+∑i=1n(∫x0x1Hyiδyi(x)dx+Hyi′δyi∣x=x1−Hyi′δyi∣x=x0−∫x0x1dHyi′dxδyi(x)dx)=H∣x=x1δx1+∑i=1n(Hyi′δyi∣x=x1−Hyi′δyi∣x=x0+∫x0x1[Hyi−dHyi′dx]δyi(x)dx)=H∣x=x1δx1+∑i=1n(Hyi′∣x=x1δyi(x1)−Hyi′∣x=x0δyi(x0)+∫x0x1[Hyi−dHyi′dx]δyi(x)dx)=H∣x=x1δx1−∑i=1nHyi′∣x=x1yi′(x1)δx1+∑i=1n(Hyi′∣x=x1(δy1)i−Hyi′∣x=x0(δy0)i+∫x0x1[Hyi−dHyi′dx]δyi(x)dx)=(H−∑i=1nHyi′yi′)∣x=x1δx1+∑i=1n(Hyi′∣x=x1(δy1)i−Hyi′∣x=x0(δy0)i+∫x0x1[Hyi−dHyi′dx]δyi(x)dx)\begin{aligned} \delta J^{\star} &= \int_{x_0}^{x_1+\delta x_1}H(x,\mathbf{y}(x)+\delta\mathbf{y}(x),\mathbf{y}'(x)+\delta\mathbf{y}'(x))\mathrm{d}x-\int_{x_0}^{x_1}H(x,\mathbf{y}(x),\mathbf{y}'(x))\mathrm{d}x \\ &= \int_{x_1}^{x_1+\delta x_1}H(x,\mathbf{y}(x)+\delta\mathbf{y}(x),\mathbf{y}'(x)+\delta\mathbf{y}'(x))\mathrm{d}x \\ &\quad +\int_{x_0}^{x_1}[H(x,\mathbf{y}(x)+\delta\mathbf{y}(x),\mathbf{y}'(x)+\delta\mathbf{y}'(x))-H(x,\mathbf{y}(x),\mathbf{y}'(x))]\mathrm{d}x \\ &= H|_{x=x_1}\delta x_1 +\sum_{i=1}^{n}\left(\int_{x_0}^{x_1}[H_{y_i}\delta y_i(x)+H_{y'_i}\delta y'_i(x)]\mathrm{d}x\right) \\ &= H|_{x=x_1}\delta x_1 +\sum_{i=1}^{n}\left(\int_{x_0}^{x_1}H_{y_i}\delta y_i(x)\mathrm{d}x +\int_{x_0}^{x_1}H_{y'_i}\delta y'_i(x)\mathrm{d}x\right) \\ &= H|_{x=x_1}\delta x_1 +\sum_{i=1}^{n}\left(\int_{x_0}^{x_1}H_{y_i}\delta y_i(x)\mathrm{d}x +\int_{x_0}^{x_1}H_{y'_i}\frac{\mathrm{d}\delta y_i(x)}{\mathrm{d}x}\mathrm{d}x\right) \\ &= H|_{x=x_1}\delta x_1 +\sum_{i=1}^{n}\left(\int_{x_0}^{x_1}H_{y_i}\delta y_i(x)\mathrm{d}x +H_{y'_i}\delta y_i|_{x=x_1}-H_{y'_i}\delta y_i|_{x=x_0}-\int_{x_0}^{x_1}\frac{\mathrm{d}H_{y'_i}}{\mathrm{d}x}\delta y_i(x)\mathrm{d}x\right) \\ &= H|_{x=x_1}\delta x_1 +\sum_{i=1}^{n}\left(H_{y'_i}\delta y_i|_{x=x_1}-H_{y'_i}\delta y_i|_{x=x_0}+\int_{x_0}^{x_1}[H_{y_i}-\frac{\mathrm{d}H_{y'_i}}{\mathrm{d}x}]\delta y_i(x)\mathrm{d}x\right) \\ &= H|_{x=x_1}\delta x_1 +\sum_{i=1}^{n}\left(H_{y'_i}|_{x=x_1}\delta y_{i}(x_1)-H_{y'_i}|_{x=x_0}\delta y_{i}(x_0)+\int_{x_0}^{x_1}[H_{y_i}-\frac{\mathrm{d}H_{y'_i}}{\mathrm{d}x}]\delta y_i(x)\mathrm{d}x\right) \\ &= H|_{x=x_1}\delta x_1 -\sum_{i=1}^{n}H_{y'_i}|_{x=x_1}y'_i(x_1)\delta x_1 \\ &\quad +\sum_{i=1}^{n}\left(H_{y'_i}|_{x=x_1}(\delta\mathbf{y}_1)_i-H_{y'_i}|_{x=x_0}(\delta\mathbf{y}_0)_i+\int_{x_0}^{x_1}[H_{y_i}-\frac{\mathrm{d}H_{y'_i}}{\mathrm{d}x}]\delta y_i(x)\mathrm{d}x\right) \\ &= \left(H -\sum_{i=1}^{n}H_{y'_i}y'_i\right)|_{x=x_1}\delta x_1 \\ &\quad +\sum_{i=1}^{n}\left(H_{y'_i}|_{x=x_1}(\delta\mathbf{y}_1)_i-H_{y'_i}|_{x=x_0}(\delta\mathbf{y}_0)_i+\int_{x_0}^{x_1}[H_{y_i}-\frac{\mathrm{d}H_{y'_i}}{\mathrm{d}x}]\delta y_i(x)\mathrm{d}x\right) \\ \end{aligned}
由于 δx1,δy(x)\delta x_1,\delta\mathbf{y}(x) 都是任意取的，互相独立，为了保证 δJ⋆=0\delta J^{\star}=0，则必有：
{(H−∑i=1nHyi′yi′)∣x=x1=0,⇐¬(δx1=0)Hyi′∣x=x0=0,⇐¬((δy0)i=0)Hyi′∣x=x1=0,⇐¬((δy1)i=0)Hyi−ddxHyi′=0\begin{cases} \left(H -\sum_{i=1}^{n}H_{y'_i}y'_i\right)|_{x=x_1}=0, & \Leftarrow \neg (\delta x_1= 0) \\ H_{y'_i}|_{x=x_0}=0,& \Leftarrow \neg ((\delta\mathbf{y}_0)_i= 0) \\ H_{y'_i}|_{x=x_1}=0,& \Leftarrow \neg((\delta\mathbf{y}_1)_i= 0) \\ H_{y_i}-\frac{\mathrm{d}}{\mathrm{d}x}H_{y'_i}=0 \end{cases}
由于起点终点位置固定，有 δy1=δy0=0\delta\mathbf{y}_1=\delta\mathbf{y}_0=\mathbf{0}，且 δx1\delta x_1 可能不是 00，故：
{(H−∑i=1nHyi′yi′)∣x=x1=0Hyi−ddxHyi′=0\begin{cases} \left(H -\sum_{i=1}^{n}H_{y'_i}y'_i\right)|_{x=x_1}=0 \\ H_{y_i}-\frac{\mathrm{d}}{\mathrm{d}x}H_{y'_i}=0 \end{cases}
由此，我们得到了欧拉-拉格朗日方程与其边界条件，代入 H=L+∑i=1mλi(x)giH=L+\sum_{i=1}^{m}\lambda_i(x)g_i 与 gi=0g_i=0，可得
(H−∑i=1nHyi′yi′)∣x=x1=(L+∑j=1mλj(x)gj−∑i=1n(L+∑j=1mλj(x)gj)yi′yi′)∣x=x1=(L−∑i=1n(Lyi′+∑j=1mλj(x)(gj)yi′)yi′)∣x=x1=(L−∑i=1nLyi′yi′−∑j=1m(λj(x)∑i=1n(gj)yi′yi′))∣x=x1\begin{aligned} &\quad\left(H -\sum_{i=1}^{n}H_{y'_i}y'_i\right)|_{x=x_1} \\ &= \left(L+\sum_{j=1}^{m}\lambda_j(x)g_j -\sum_{i=1}^{n}(L+\sum_{j=1}^{m}\lambda_j(x)g_j)_{y'_i}y'_i\right)|_{x=x_1} \\ &= \left(L -\sum_{i=1}^{n}(L_{y'_i}+\sum_{j=1}^{m}\lambda_j(x)(g_j)_{y'_i})y'_i\right)|_{x=x_1} \\ &= \left(L -\sum_{i=1}^{n}L_{y'_i}y'_i-\sum_{j=1}^{m}\left(\lambda_j(x)\sum_{i=1}^{n}(g_j)_{y'_i}y'_i\right)\right)|_{x=x_1} \\ \end{aligned} Hyi−ddxHyi′=(L+∑j=1mλj(x)gj)yi−ddx(L+∑j=1mλj(x)gj)yi′=Lyi+∑j=1mλj(x)(gj)yi−ddxLyi′−∑j=1mddx(λj(x)(gj)yi′)=Lyi−ddxLyi′+∑j=1m(λj(x)(gj)yi−ddx(λj(x)(gj)yi′))\begin{aligned} &\quad H_{y_i}-\frac{\mathrm{d}}{\mathrm{d}x}H_{y'_i} \\ &= (L+\sum_{j=1}^{m}\lambda_j(x)g_j)_{y_i}-\frac{\mathrm{d}}{\mathrm{d}x}(L+\sum_{j=1}^{m}\lambda_j(x)g_j)_{y'_i} \\ &= L_{y_i}+\sum_{j=1}^{m}\lambda_j(x)(g_j)_{y_i}-\frac{\mathrm{d}}{\mathrm{d}x}L_{y'_i}-\sum_{j=1}^{m}\frac{\mathrm{d}}{\mathrm{d}x}(\lambda_j(x)(g_j)_{y'_i}) \\ &= L_{y_i}-\frac{\mathrm{d}}{\mathrm{d}x}L_{y'_i}+\sum_{j=1}^{m}\left(\lambda_j(x)(g_j)_{y_i}-\frac{\mathrm{d}}{\mathrm{d}x}(\lambda_j(x)(g_j)_{y'_i})\right) \\ \end{aligned} {(L−∑i=1nLyi′yi′−∑j=1m(λj(x)∑i=1n(gj)yi′yi′))∣x=x1=0Lyi−ddxLyi′+∑j=1m(λj(x)(gj)yi−ddx(λj(x)(gj)yi′))=0gi=0\begin{cases} \left(L -\sum_{i=1}^{n}L_{y'_i}y'_i-\sum_{j=1}^{m}\left(\lambda_j(x)\sum_{i=1}^{n}(g_j)_{y'_i}y'_i\right)\right)|_{x=x_1}=0 \\ L_{y_i}-\frac{\mathrm{d}}{\mathrm{d}x}L_{y'_i}+\sum_{j=1}^{m}\left(\lambda_j(x)(g_j)_{y_i}-\frac{\mathrm{d}}{\mathrm{d}x}(\lambda_j(x)(g_j)_{y'_i})\right)=0 \\ g_i=0 \end{cases}
这样我们就得到了有约束情况下的欧拉-拉格朗日方程。
从此出发，我们从费马原理求光线方程，对于费马原理：
δS=δ∫ABnds=0\delta S=\delta \int_{A}^B n\mathrm{d}s=0
我们可以知道 L=n(r)L=n(\mathbf{r}) 与 r′\mathbf{r}' 无关，设所处空间维度数为 kk，r(s)=(x1(s),x2(s),…xk(s))\mathbf{r}(s)=(x_1(s),x_2(s),\ldots x_k(s)) 为位矢，由于曲线是以弧长 ss 为自变量，所以有一个约束：
g(r′(s))=∑i=1k(dxids)2−1=∑i=1kx′i2−1=(drds)2−1=0g(\mathbf{r}'(s))=\sum_{i=1}^{k}\left(\frac{\mathrm{d}x_i}{\mathrm{d}s}\right)^2-1=\sum_{i=1}^{k} {x'}_i^2-1=\left(\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}s}\right)^2-1=0
计算 dgds=0\frac{\mathrm{d}g}{\mathrm{d}s}=0 我们还可以得到
dgds=2∑i=1kdxidsd2xids2=0\frac{\mathrm{d}g}{\mathrm{d}s}=2\sum_{i=1}^{k}\frac{\mathrm{d}x_i}{\mathrm{d}s}\frac{\mathrm{d}^2 x_i}{\mathrm{d}s^2}=0
那么对应的欧拉-拉格朗日方程为：
(n−∑i=1knxi′xi′−(λ(s)∑i=1kgxi′xi′))∣s=s1=(n−(λ(s)∑i=1k(2xi′)xi′))∣s=s1=(n−λ(s)(2+2g))∣s=s1=(n−2λ(s))∣s=s1\begin{aligned} &\quad\left(n -\sum_{i=1}^{k}n_{x'_i}x'_i-\left(\lambda(s)\sum_{i=1}^{k}g_{x'_i}x'_i\right)\right)|_{s=s_1} \\ &= \left(n -\left(\lambda(s)\sum_{i=1}^{k}(2x'_i)x'_i\right)\right)|_{s=s_1} \\ &= \left(n -\lambda(s)(2+2g)\right)|_{s=s_1} \\ &= \left(n -2\lambda(s)\right)|_{s=s_1} \\ \end{aligned} nxi−ddsnxi′+(λ(s)gxi−dds(λ(s)gxi′))=nxi−dds0+(λ(s)0−dds(λ(s)(2xi′)))=∂n∂xi−dds(2λ(s)xi′)\begin{aligned} &\quad n_{x_i}-\frac{\mathrm{d}}{\mathrm{d}s}n_{x'_i}+\left(\lambda(s)g_{x_i}-\frac{\mathrm{d}}{\mathrm{d}s}(\lambda(s)g_{x'_i})\right) \\ &= n_{x_i}-\frac{\mathrm{d}}{\mathrm{d}s}0+\left(\lambda(s)0-\frac{\mathrm{d}}{\mathrm{d}s}(\lambda(s)(2x'_i))\right) \\ &= \frac{\partial n}{\partial x_i}-\frac{\mathrm{d}}{\mathrm{d}s}(2\lambda(s)x'_i) \\ \end{aligned} {(n−2λ(s))∣s=s1=0∂n∂xi−dds(2λ(s)xi′)=0g=∑i=1kx′i2−1=0\begin{cases} \left(n -2\lambda(s)\right)|_{s=s_1}=0 \\ \frac{\partial n}{\partial x_i}-\frac{\mathrm{d}}{\mathrm{d}s}(2\lambda(s)x'_i)=0 \\ g=\sum_{i=1}^{k} {x'}_i^2-1=0 \end{cases}
我们可以把中间的式子写作矢量式，通过给第 ii 个方程乘以 xix_i 方向的单位向量 x^i\hat{x}_i 再加和，有：
0=∑i=1k(∂n∂xi−dds(2λ(s)xi′))x^i=∑i=1k∂n∂xix^i−dds(2λ(s)d∑i=1kxix^ids)=∇n−dds(2λ(s)dr(s)ds)\begin{aligned} \mathbf{0} &=\sum_{i=1}^{k}\left(\frac{\partial n}{\partial x_i}-\frac{\mathrm{d}}{\mathrm{d}s}(2\lambda(s)x'_i)\right)\hat{x}_i \\ &=\sum_{i=1}^{k}\frac{\partial n}{\partial x_i}\hat{x}_i-\frac{\mathrm{d}}{\mathrm{d}s}(2\lambda(s)\frac{\mathrm{d}\sum_{i=1}^{k}x_i\hat{x}_i}{\mathrm{d}s}) \\ &=\nabla n-\frac{\mathrm{d}}{\mathrm{d}s}(2\lambda(s)\frac{\mathrm{d}\mathbf{r}(s)}{\mathrm{d}s}) \\ \end{aligned}
我们考虑沿着曲线 r(s)\mathbf{r}(s) 上 λ(s)\lambda(s) 的变化，通过把得到的向量方程向曲线切向投影，即把中间的式子中第 ii 个方程乘以 dxids\frac{\mathrm{d}x_i}{\mathrm{d}s} 再加和：
0=∑i=1k(∂n∂xi−dds(2λ(s)xi′))dxids=∑i=1k∂n∂xidxids−∑i=1kdds(2λ(s)dxids)dxids=dnds−d(2λ(s))ds∑i=1k(dxids)2−2λ(s)∑i=1kd2xids2dxids=dnds−d(2λ(s))ds\begin{aligned} 0 &=\sum_{i=1}^{k}\left(\frac{\partial n}{\partial x_i}-\frac{\mathrm{d}}{\mathrm{d}s}(2\lambda(s)x'_i)\right)\frac{\mathrm{d}x_i}{\mathrm{d}s} \\ &=\sum_{i=1}^{k}\frac{\partial n}{\partial x_i}\frac{\mathrm{d}x_i}{\mathrm{d}s}-\sum_{i=1}^{k}\frac{\mathrm{d}}{\mathrm{d}s}(2\lambda(s)\frac{\mathrm{d}x_i}{\mathrm{d}s})\frac{\mathrm{d}x_i}{\mathrm{d}s} \\ &=\frac{\mathrm{d}n}{\mathrm{d}s}-\frac{\mathrm{d}(2\lambda(s))}{\mathrm{d}s}\sum_{i=1}^{k}\left(\frac{\mathrm{d}x_i}{\mathrm{d}s}\right)^2-2\lambda(s)\sum_{i=1}^{k}\frac{\mathrm{d}^2 x_i}{\mathrm{d}s^2}\frac{\mathrm{d}x_i}{\mathrm{d}s} \\ &=\frac{\mathrm{d}n}{\mathrm{d}s}-\frac{\mathrm{d}(2\lambda(s))}{\mathrm{d}s}\\ \end{aligned}
解上述方程，我们可以知道 2λ(s)=n(r(s))+C2\lambda(s)=n(\mathbf{r}(s))+C，CC 是某个常数。根据边界条件 (n−2λ(s))∣s=s1=0\left(n -2\lambda(s)\right)|_{s=s_1}=0 可以知道常数 C=0C=0，于是 2λ(s)=n(r(s))2\lambda(s)=n(\mathbf{r}(s)) ，代入矢量式，可得：
dds(ndr(s)ds)=∇n\frac{\mathrm{d}}{\mathrm{d}s}(n\frac{\mathrm{d}\mathbf{r}(s)}{\mathrm{d}s})=\nabla n
于是我们得到了光线方程。

(2) 从光线方程到以 xx 为自变量的微分方程

我们从光线方程出发，考虑空间维度 k=2,x1=x,x2=yk=2,x_1=x,x_2=y，其在 xx 方向投影（不考虑 yy 方向上是因为最终得到的方程是相同的），有：
dds(ndxds)=∂n∂x\frac{\mathrm{d}}{\mathrm{d}s}(n\frac{\mathrm{d}x}{\mathrm{d}s})=\frac{\partial n}{\partial x}
代入弧长 ds=1+(dydx)2dx\mathrm{d}s=\sqrt{1+\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2} \mathrm{d}x 我们可以得到：
11+(dydx)2ddx(n11+(dydx)2)=∂n∂x1(1+(dydx)2)12−ndydxd2ydx2+dndx+dndx(dydx)2(1+(dydx)2)32=∂n∂x−ndydxd2ydx21+(dydx)2+dndx=∂n∂x(1+(dydx)2)−ndydxd2ydx21+(dydx)2+∂n∂x+∂n∂ydydx=∂n∂x(1+(dydx)2)−∂n∂x(dydx)2+∂n∂ydydx=ndydxd2ydx21+(dydx)2∂n∂y−∂n∂xdydx=nd2ydx21+(dydx)21+(dydx)2n(∂n∂y−∂n∂xdydx)=d2ydx2y′′=1+y′2n(∂n∂y−y′∂n∂x)\begin{aligned} \frac{1}{\sqrt{1+\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2}}\frac{\mathrm{d}}{\mathrm{d}x}(n\frac{1}{\sqrt{1+\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2}})&=\frac{\partial n}{\partial x} \\ \frac{1}{\left(1+\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2\right)^{\frac{1}{2}}}\frac{-n\frac{\mathrm{d}y}{\mathrm{d}x}\frac{\mathrm{d}^2 y}{\mathrm{d}x^2}+\frac{\mathrm{d}n}{\mathrm{d}x}+\frac{\mathrm{d}n}{\mathrm{d}x}\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2}{\left(1+\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2\right)^{\frac{3}{2}}}&=\frac{\partial n}{\partial x} \\ \frac{-n\frac{\mathrm{d}y}{\mathrm{d}x}\frac{\mathrm{d}^2 y}{\mathrm{d}x^2}}{1+\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2}+\frac{\mathrm{d}n}{\mathrm{d}x}&=\frac{\partial n}{\partial x}\left(1+\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2\right) \\ \frac{-n\frac{\mathrm{d}y}{\mathrm{d}x}\frac{\mathrm{d}^2 y}{\mathrm{d}x^2}}{1+\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2}+\frac{\partial n}{\partial x}+\frac{\partial n}{\partial y}\frac{\mathrm{d}y}{\mathrm{d}x}&=\frac{\partial n}{\partial x}\left(1+\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2\right) \\ -\frac{\partial n}{\partial x}\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2+\frac{\partial n}{\partial y}\frac{\mathrm{d}y}{\mathrm{d}x}&=\frac{n\frac{\mathrm{d}y}{\mathrm{d}x}\frac{\mathrm{d}^2 y}{\mathrm{d}x^2}}{1+\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2} \\ \frac{\partial n}{\partial y}-\frac{\partial n}{\partial x}\frac{\mathrm{d}y}{\mathrm{d}x}&=\frac{n\frac{\mathrm{d}^2 y}{\mathrm{d}x^2}}{1+\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2} \\ \frac{1+\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2}{n}\left(\frac{\partial n}{\partial y}-\frac{\partial n}{\partial x}\frac{\mathrm{d}y}{\mathrm{d}x}\right)&=\frac{\mathrm{d}^2 y}{\mathrm{d}x^2} \\ y''&=\frac{1+y'^2}{n}\left(\frac{\partial n}{\partial y}-y'\frac{\partial n}{\partial x}\right) \\ \end{aligned}
于是得到以 xx 为自变量的微分方程。
如果要以 ss 为自变量，则有：
dds(ndxids)=∂n∂xidndsdxids+ndxi2ds2=∂n∂xi(∑j=1k∂n∂xjdxjds)dxids+ndxi2ds2=∂n∂xidxi2ds2=1n(∂n∂xi−(∑j=1k∂n∂xjdxjds)dxids)\begin{aligned} \frac{\mathrm{d}}{\mathrm{d}s}(n\frac{\mathrm{d}x_i}{\mathrm{d}s})&=\frac{\partial n}{\partial x_i} \\ \frac{\mathrm{d}n}{\mathrm{d}s}\frac{\mathrm{d}x_i}{\mathrm{d}s}+n\frac{\mathrm{d}x_i^2}{\mathrm{d}s^2}&=\frac{\partial n}{\partial x_i} \\ \left(\sum_{j=1}^k\frac{\partial n}{\partial x_j}\frac{\mathrm{d}x_j}{\mathrm{d}s}\right)\frac{\mathrm{d}x_i}{\mathrm{d}s}+n\frac{\mathrm{d}x_i^2}{\mathrm{d}s^2}&=\frac{\partial n}{\partial x_i} \\ \frac{\mathrm{d}x_i^2}{\mathrm{d}s^2}&=\frac{1}{n}\left(\frac{\partial n}{\partial x_i}-\left(\sum_{j=1}^k\frac{\partial n}{\partial x_j}\frac{\mathrm{d}x_j}{\mathrm{d}s}\right)\frac{\mathrm{d}x_i}{\mathrm{d}s}\right) \\ \end{aligned}
解上述微分方程时，对于约束成立的初始条件，约束 ∑i=1k(dxids)2−1=0\sum_{i=1}^{k}\left(\frac{\mathrm{d}x_i}{\mathrm{d}s}\right)^2-1=0 应自动保持成立，如存在误差可在每一步对一阶导向量 drds\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}s} 进行归一化。

(3) 从折射定律进行微元分析得到光子角动量守恒

不妨设光线随极坐标的 ϕ\phi 变大向外行进(由于光路可逆和对称性可知其他情况也成立)，即 ri+1>rir_{i+1}>r_i，将介质分为无穷薄的多层球壳，每层介质的折射率为 nin_i，其下界面到极坐标的极点距离为 rir_i，光线从下界面的出射角为 θi\theta_i，光线与下界面交点到极坐标的极点连线与极坐标的极轴夹角 ϕi\phi_i，那么 dϕ=ϕi+1−ϕi>0\mathrm{d}\phi=\phi_{i+1}-\phi_{i}>0 也是无穷小，有折射定律（计算中只保留到一阶无穷小）：
ni+1sin⁡θn+1=nisin⁡(θn−dϕ)ri+1ni+1sin⁡θn+1=ri+1nisin⁡(θn−dϕ)ri+1ni+1sin⁡θn+1=ni(ri+1sin⁡θn−cos⁡θnri+1dϕ)ri+1ni+1sin⁡θn+1=nisin⁡θn(ri+1−cot⁡θnri+1dϕ)ri+1ni+1sin⁡θn+1=nisin⁡θn(ri+1−(ri+1−ri))ri+1ni+1sin⁡θn+1=rinisin⁡θn\begin{aligned} n_{i+1}\sin\theta_{n+1}&=n_{i}\sin(\theta_{n}-\mathrm{d}\phi) \\ r_{i+1}n_{i+1}\sin\theta_{n+1}&=r_{i+1}n_{i}\sin(\theta_{n}-\mathrm{d}\phi) \\ r_{i+1}n_{i+1}\sin\theta_{n+1}&=n_{i}(r_{i+1}\sin\theta_{n}-\cos\theta_{n}r_{i+1}\mathrm{d}\phi) \\ r_{i+1}n_{i+1}\sin\theta_{n+1}&=n_{i}\sin\theta_{n}(r_{i+1}-\cot\theta_{n}r_{i+1}\mathrm{d}\phi) \\ r_{i+1}n_{i+1}\sin\theta_{n+1}&=n_{i}\sin\theta_{n}(r_{i+1}-(r_{i+1}-r_i)) \\ r_{i+1}n_{i+1}\sin\theta_{n+1}&=r_in_{i}\sin\theta_{n} \\ \end{aligned}
于是可知
rn(r)sin⁡θ=constrn(r)\sin{\theta}=const
实际上介质中光子动量为 p=nℏωcp=n\frac{\hbar\omega}{c}，可以看到是正比于 nn 的，那么上式可以写成
rpsin⁡θ=constrp\sin{\theta}=const
在二维极坐标中，rpsin⁡θ=∣r×p∣=∣L∣=constrp\sin{\theta}=|\mathbf{r}\times \mathbf{p}|=|\mathbf{L}|=const，L\mathbf{L} 为角动量，光子角动量守恒。

最后附上代码

CPP
// Radioactive Islands.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
#include <iostream>
#include <vector>
#include <algorithm>
#include <iomanip>
#define _USE_MATH_DEFINES
#include <math.h>
#include <float.h>
using namespace std;
//Bogacki和Shampine的5(4)格式,在mathematica中可使用NDSolve`EmbeddedExplicitRungeKuttaCoefficients[5, Infinity]得到
constexpr int RKp = 5;
constexpr int RKn = 8;
constexpr double RKa[] = { 0.16666666666666666667,0.074074074074074074074,0.14814814814814814815,0.13338192419825072886,-0.47230320699708454810,0.76749271137026239067,0.22895622895622895623,-0.36363636363636363636,0.29370629370629370629,0.50764050764050764051,0.026500355113636363636,0.23011363636363636364,0.10772747760052447552,0.16021907779720279720,0.22543945312500000000,0.18159821692916505081,-0.38707982536247294745,0.41288268560074054269,0.64099131912534967441,-1.0121439383514517683,1.1637515420586694478,0.072792658730158730159,0,0.28662437773692472941,0.19513621794871794872,0.0086383928571428571429,0.35956558061821219716,0.077242772108843537415 };
constexpr double RKc[] = { 0,0.16666666666666666667,0.22222222222222222222,0.42857142857142857143,0.66666666666666666667,0.75000000000000000000,1.0000000000000000000,1.0000000000000000000 };
constexpr double RKbErr[] = { 0.0019478942125547063819,0,-0.0090486991861521936710,0.015478743126634073136,-0.021222718253968253968,0.013276905111429694492,0.0054803707701381459657,-0.0059124957806361723368 };
//常微分方程的状态，视作向量，分别是已受到辐射量，y一阶导，y
struct State
{
    double score, dy, y;
    public:
        State()
        {
            this->score = 0;
            this->dy = 0;
            this->y = 0;
        }
        State operator+(const State& b)
        {
            State state;
            state.score = this->score + b.score;
            state.dy = this->dy + b.dy;
            state.y = this->y + b.y;
            return state;
        }
        State operator*(double b)
        {
            State state;
            state.score = this->score * b;
            state.dy = this->dy * b;
            state.y = this->y * b;
            return state;
        }
        double norm()
        {
            return sqrt(score* score+dy * dy + y * y);
        }
        void operator+=(const State& b)
        {
            this->score+= b.score;
            this->dy+= b.dy;
            this->y+= b.y;
        }
        void operator*=(double b)
        {
            this->score *= b;
            this->dy *= b;
            this->y *= b;
        }
};
//求解轨迹
class Trajectory
{
public:
    Trajectory(double A, double B, vector<double> C);
    ~Trajectory();
    double solve(double theta);
    void get_y0_int();
    State y;
    double Tmax,y0;
    int y0_int;
//private:
    double tol = 0.0000001,min_turn=0.001;
    double A, B;
    int counter;
    vector<double> C;
    State f(State y,double x);
    bool dy_divergeQ(State y, double x, double h);
    bool fall_into_pointQ(State y, double x);
    bool interval_sol_exist(int interval);
};
//如果两点挨得足够近，则通过其中的必然不是最优解
bool Trajectory::interval_sol_exist(int interval)
{
    if (interval==0 || interval== C.size())
    {
        return true;
    }
    double minT{ sqrt(400.0 + (A - B) * (A - B)) };
    for (int i = 0; i < C.size(); i++)
    {
        double d = max(abs(C[interval] - C[i]), abs(C[interval - 1] - C[i]));
        minT += max(2.0 / d - 0.1 - 0.1, 0.0);
    }
    return minT < Tmax;
}
//微分方程右端的f函数
State Trajectory::f(State y,double x)
{
    State dy;
    dy.y = y.dy;
    dy.dy = 0.0;
    dy.score = 1.0;
    for (int i = 0; i < C.size(); i++)
    {
        double temp = 1.0 / (x * x + (y.y - C[i]) * (y.y - C[i]));
        dy.dy += (-(y.y - C[i]) + y.dy * x) * temp*temp;
        dy.score += temp;
    }
    dy.dy *= (1 + y.dy * y.dy) * 2.0/ dy.score;
    dy.score *= sqrt(1 + y.dy * y.dy);
    return dy;
}
//将y(0)转换成处于哪个区间内
void Trajectory::get_y0_int()
{
    if (y0<=C[0])
    {
        y0_int = 0;
        return;
    }
    if (C.back() < y0)
    {
        y0_int = C.size();
        return;
    }
    int min{ 1 }, max{ int(C.size()) - 1 };
    while (true)
    {
        int mid = (min + max) / 2;
        if (y0 <= C[mid-1])
        {
            max = mid - 1;
            continue;
        }
        if (C[mid] < y0)
        {
            min = mid + 1;
            continue;
        }
        y0_int = mid;
        return;
    }
}
//检查轨迹是否落入某一点
bool Trajectory::fall_into_pointQ(State y, double x)
{
    double n = 1.0;
    double maxn{ 0 };
    int maxi{ -1 };
    for (int i = 0; i < C.size(); i++)
    {
        double temp = 1.0 / (x * x + (y.y - C[i]) * (y.y - C[i]));
        if (maxn<temp)
        {
            maxn = temp;
            maxi = i;
        }
        n += temp;
    }
    if((n-maxn)/n<0.01)
    {
        double d = sqrt(x * x + (y.y - C[maxi]) * (y.y - C[maxi]));
        double d0 = 1.0 / d;
        double d0x{ -x * d0 }, d0y{ -(y.y - C[maxi]) * d0 };
        double dy0 = 1.0 / sqrt(1 + y.dy * y.dy);
        double dyx{ dy0 }, dyy{ y.dy*dy0 };
        double cos0 = d0x * dyx + d0y * dyy;
        if (cos0>0.0 && 1.0 / (n * d * sqrt(1 - cos0 * cos0)) > 2.0 * d)
        {
            y0 = C[maxi];
            y0_int = maxi;
            return true;
        }
    }
    return false;
}
//检查是否y导数发散到无穷，即轨迹竖直
bool Trajectory::dy_divergeQ(State y, double x, double h)
{
    double theta = atan(1.0 / y.dy);
    double sintheta{ 1.0 / sqrt(1 + y.dy * y.dy) }, costheta{ y.dy / sqrt(1 + y.dy * y.dy) };
    double R = 0;
    double n = 1;
    for (int i = 0; i < C.size(); i++)
    {
        double temp = 1.0 / (x * x + (y.y - C[i]) * (y.y - C[i]));
        R += (-(y.y - C[i]) * sintheta + x * costheta) * temp * temp * 2.0;
        n += temp;
    }
    R /= n;
    R = 1.0 / R;
    return y.dy * R > 0 && theta * abs(R) < min_turn && abs(R) * (1 - costheta) < h * 0.5;
}
constexpr double xmin = -10;
constexpr double xmax = 10;
//求解轨迹
double Trajectory::solve(double theta)
{
    vector<State> k(RKn);
    double h = (xmax - xmin) * tol;
    y.y = A;
    y.dy = tan(theta);
    y.score = 0;
    double x = xmin;
    bool sol_notfound = true;
    bool zero_notreached = true;
    bool zero_reached = false;
    bool zero_flag = false;
    k[0] = f(y, x);
    counter = 0;
    while (sol_notfound)
    {
        if (zero_notreached && x + h >= 0)//x=0
        {
            zero_notreached = false;
            zero_flag = true;
            h = - x;
        }
        if (x+h>=xmax)//x=xmax
        {
            sol_notfound = false;
            h = xmax - x;
        }
        if (abs(y.y)-10 > 0.5*Tmax) //y出界
        {
            if (zero_notreached)
            {
                y0 = y.y;
                get_y0_int();
            }
            if (y.y > 0)
            {
                return DBL_MAX;
            }
            else
            {
                return -DBL_MAX;
            }

        }
        if (zero_notreached && fall_into_pointQ(y, x))//落入某一点
        {
            return DBL_MAX;

        }
        if (zero_reached && y.score + sqrt((y.y - B) * (y.y - B) + (x - xmax) * (x - xmax)) > Tmax)//过0且过长
        {
            if (y.dy > 0)
            {
                return DBL_MAX;
            }
            else
            {
                return -DBL_MAX;
            }

        }
        if (dy_divergeQ(y,x,xmax-x))//y方向发散
        {
            if (zero_notreached)
            {
                y0 = y.y;
                get_y0_int();
            }
            if (y.dy>0)
            {
                return DBL_MAX;
            }
            else
            {
                return -DBL_MAX;
            }
            
        }
        //RK迭代
        int RKa_offset = 0;
        for (int i = 1; i < RKn -1; i++)
        {
            State g;
            for (int j = 0; j < i; j++)
            {
                g += k[j] * RKa[j + RKa_offset];
            }
            g *= h;
            g += y;
            k[i] = f(g, x + RKc[i]*h);
            RKa_offset += i;
        }
        //First Same As Last （FSAL）技巧
        State y_next;
        for (int j = 0; j < RKn - 1; j++)
        {
            y_next += k[j] * RKa[j + RKa_offset];
        }
        y_next *= h;
        y_next += y;
        k[RKn - 1] = f(y_next, x + h);
        //误差估计
        State err;
        for (int i = 0; i < RKn; i++)
        {
            err += k[i] * RKbErr[i];
        }
        err *= h;
        k[0] = k[RKn - 1];
        y = y_next;
        x += h;
        //调整步长
        h *= pow((tol / err.norm()), 1.0 / double(RKp));
        if (zero_flag)
        {
            y0 = y.y;
            get_y0_int();
            zero_flag = false;
            zero_reached = true;
        }
        counter++;
    }
    if (y.y>0.5*Tmax)
    {
        return DBL_MAX;
    }
    if (y.y < -0.5 * Tmax)
    {
        return -DBL_MAX;
    }
    return y.y - B;
}

Trajectory::Trajectory(double A, double B, vector<double>C)
{
    this->A = A;
    this->B = B;
    this->C = C;
    sort(this->C.begin(), this->C.end());
    Tmax = (min(20 - A - B, 20 + A + B) + M_PI * 10) * (1 + 0.01 * C.size());
}

Trajectory::~Trajectory()
{

}
//在某区间中的一点
struct Point
{
    double x, fx;
    int interval;
public:
    Point(double x,double fx,int interval)
    {
        this->x = x;
        this->fx = fx;
        this->interval = interval;
    }
    friend void swap(Point& first, Point& second) noexcept {
        swap(first.x, second.x);
        swap(first.fx, second.fx);
        swap(first.interval, second.interval);
    }
    bool operator >(const Point& d)
    {
        return x > d.x;
    }
};
//只用于Brent法
struct Point2d
{
    double x, y;
public:
    Point2d(double x, double y)
    {
        this->x = x;
        this->y = y;
    }
    friend void swap(Point2d& first, Point2d& second) noexcept {
        swap(first.x, second.x);
        swap(first.y, second.y);
    }
};
//Brent法解方程
double findRoot_Brent(pair<Point2d, Point2d> interval, double epsif,double epsix,double delta,int itermax,Trajectory& tr)
{
    Point2d a(interval.first.x, interval.first.y), b(interval.second.x, interval.second.y);
    if (a.y==0)
    {
        return a.x;
    }
    if (b.y == 0)
    {
        return b.x;
    }
    if (a.y*b.y >0)
    {
        return DBL_MAX;
    }
    if (abs(a.y)<abs(b.y))
    {
        swap(a, b);
    }
    Point2d c(a);
    bool mflag = true;
    Point2d s(0, 0);
    double d = 0;
    int iter = 0;
    while ((abs(b.y)>epsif || abs(b.x-a.x)>epsix )&&(iter<=itermax))
    {
        iter++;
        if ((! c.y==a.y)&& (!c.y == b.y))
        {
            s.x = a.x * b.y * c.y / ((a.y - b.y) * (a.y - c.y)) + b.x * a.y * c.y / ((b.y - a.y) * (b.y - c.y)) + c.x * a.y * b.y / ((c.y - a.y) * (c.y - b.y));
        }
        else
        {
            s.x = b.x - b.y * (b.x - a.x) / (b.y - a.y);
        }
        if (((s.x - b.x) * (s.x - (3.0 * a.x + b.x) * 0.25) > 0) ||
            (mflag && abs(s.x - b.x) >= 0.5 * abs(b.x - c.x)) ||
            ((!mflag) && abs(s.x - b.x) >= 0.5 * abs(c.x - d)) ||
            (mflag && abs(b.x - c.x) < delta) ||
            ((!mflag) && abs(c.x - d) < delta))
        {
            s.x = 0.5 * (a.x + b.x);
            mflag = true;
        }
        else
        {
            mflag = false;
        }
        s.y = tr.solve(s.x);
        d = c.x;
        c = b;
        if (s.y*a.y<0)
        {
            b = s;
        }
        else
        {
            a = s;
        }
        if (abs(a.y) < abs(b.y))
        {
            swap(a, b);
        }
    }
    return b.x;
}

constexpr double thetaTol=0.00000001;
//求出最低辐射量
double leastScore(double A, double B, vector<double> C)
{
    Trajectory tr(A, B, C);
    vector<pair<Point, Point>> interval(tr.C.size()+1, pair<Point, Point>(Point(-M_PI_2, -DBL_MAX,0), Point(M_PI_2, DBL_MAX, int(tr.C.size()))));
    vector<bool> interval_sol_exist(tr.C.size() + 1, true);
    //排除绝对无解的区间
    for (int i = 0; i < interval_sol_exist.size(); i++)
    {
        interval_sol_exist[i] = tr.interval_sol_exist(i);
    }
    //对每个区间逐次进行二分法使得interval[i]里存的俩点处于第i个区间内
    for (int i = 0; i < interval.size(); i++)
    {
        if (!interval_sol_exist[i])
        {
            continue;
        }
        while ((interval[i].first.interval != i || interval[i].second.interval!=i) && (interval[i].second.x - interval[i].first.x)> thetaTol)
        {
            double mid = (interval[i].first.x + interval[i].second.x) * 0.5;
            double fmid=tr.solve(mid);
            for (int j = i; j < tr.y0_int; j++)
            {
                if (interval[j].second.x>mid)
                {
                    interval[j].second.x = mid;
                    interval[j].second.fx = fmid;
                    interval[j].second.interval = tr.y0_int;
                }
            }
            if (tr.y0_int>= i)
            {
                if (fmid > 0.0)
                {
                    if (interval[tr.y0_int].second.x > mid)
                    {
                        interval[tr.y0_int].second.x = mid;
                        interval[tr.y0_int].second.fx = fmid;
                        interval[tr.y0_int].second.interval = tr.y0_int;
                    }
                }
                else
                {
                    if (interval[tr.y0_int].first.x < mid)
                    {
                        interval[tr.y0_int].first.x = mid;
                        interval[tr.y0_int].first.fx = fmid;
                        interval[tr.y0_int].first.interval = tr.y0_int;
                    }
                }
            }
            for (int j = tr.y0_int+1; j < interval.size(); j++)
            {
                if (interval[j].first.x < mid)
                {
                    interval[j].first.x = mid;
                    interval[j].first.fx = fmid;
                    interval[j].first.interval = tr.y0_int;
                }
            }
        }
        if (interval[i].first.interval != i || interval[i].second.interval != i)
        {
            interval_sol_exist[i] = false;
        }
    }
    //对第i个区间内进行二分，直到函数值tr.solve(mid)不是正负DBL_MAX(代表无解)
    for (int i = 0; i < interval.size(); i++)
    {
        if (!interval_sol_exist[i])
        {
            continue;
        }
        while (max(abs(interval[i].first.fx), abs(interval[i].second.fx)) == DBL_MAX && (interval[i].second.x - interval[i].first.x) > thetaTol)
        {
            double mid = (interval[i].first.x + interval[i].second.x) * 0.5;
            double fmid = tr.solve(mid);
            if (fmid > 0.0)
            {
                interval[i].second.x = mid;
                interval[i].second.fx = fmid;
                interval[i].second.interval = tr.y0_int;
            }
            else
            {
                interval[i].first.x = mid;
                interval[i].first.fx = fmid;
                interval[i].first.interval = tr.y0_int;
            }
        }
        if (max(abs(interval[i].first.fx), abs(interval[i].second.fx)) == DBL_MAX)
        {
            interval_sol_exist[i] = false;
        }
    }
    //对每个区间进行求根，比较得到最少辐射量
    double least{ DBL_MAX };
    for (int i = 0; i < interval.size(); i++)
    {
        if (!interval_sol_exist[i])
        {
            continue;
        }
        double root = findRoot_Brent(pair<Point2d, Point2d>(Point2d(interval[i].first.x, interval[i].first.fx), Point2d(interval[i].second.x, interval[i].second.fx)), 10, thetaTol, 0.00000001, 100, tr);
        double froot=tr.solve(root);
        least = min(least, tr.y.score);
    }
    return least;
}

int main() {
    int n;
    cin >> n;
    vector<int> N(n);
    vector<double> A(n), B(n),result(n);
    vector<vector<double>> C(n, vector<double>());
    for (int i = 0; i < n; i++)
    {
        cin >> N[i] >> A[i] >> B[i];
        for (int j = 0; j < N[i]; j++)
        {
            double temp;
            cin >> temp;
            C[i].push_back(temp);
        }
    }

    for (int i = 0; i < n; i++)
    {
        result[i] = leastScore(A[i], B[i], C[i]);
    }
    cout << setprecision(15);
    for (int i = 0; i < n; i++)
    {
        cout << "Case #" << i + 1 << ": " << result[i] << '\n';
    }
    return 0;
}

评论

共 2 条评论，欢迎与作者交流。

正在加载评论...