使用 RK 法求解变分问题所对应的微分方程：P13210 [GCJ 2016 Finals] Radioactive Islands
洛谷(P13210) Radioactive Islands ：使用 RK 法求解变分问题所对应的微分方程
相比于其他方法，此方法计算速度可以较快，使用最后给出的代码，对于此题目的数据集，总耗时约
起点终点已确定，题目要求中，要最小化的泛函是：
我们假设
其最小值必然满足变分为
可以看到，如果我们令
那么费马原理就和原问题所需最小化函数的变分为
1. 光线方程推导
首先我们规定符号：
从一般的光线方程出发（此方程也可以从费马原理开始使用欧拉-拉格朗日方程得到，要注意的是有半完整约束
其中
注：实际上以弧长
折射率的偏导也不难求：
高阶方程皆可降次，只需把
我们在计算路径的过程中，可以顺便算出所受辐射量，只需要把
最终的微分方程组为：
2. 显式 Runge-Kutta 法解微分方程
我这里用的是 Bogacki 和 Shampine 的
NDSolve`EmbeddedExplicitRungeKuttaCoefficients[5, Infinity]这是一种具有嵌入对的 RK 法，也就是计算一步可以得到一个
以此在保证误差的情况下尽量减少迭代步数。
3. 微分方程的停止条件
首先是当
另外我们需要计算
另外有几种情况也是需要停止的，首先我们对最短路径的所受辐射上界作一估计：
这相当于沿着
(1) y 超限
那么当
(2) J 超限
另外由于我们在计算路径的过程中可以计算已受辐射量，那么可以对现在路径受到辐射的下界进行一个估计：
如果上式不成立且
另外两种停止条件没那么直观：
(3) y' 发散
第一种情况，在初始角度随意的情况下，光线可能往回拐，这会导致
判断这种情况的办法是计算当前的光线与
即圆的最右侧在距离当前位置不到
(4) 落入辐射源
第二种情况，在初始角度随意的情况下，根据微分方程演化的光线有可能撞上某个辐射源所在的点，这时所受辐射量会发散到无穷，这显然不是所受辐射最小路径。
我们从极坐标的角度来了解这件事，当我们距离某个辐射源足够近的时候，
其中
进一步的，代入
负号是因为我们假设光线向内运动，在
如果方程右侧大于
4. 将问题转化为求解非线性方程
直到现在我们讨论的都是在给定初始点 A 与初始光线角度
的零点即为我们所求的使得光线正好经过 B 点的
5. g(\theta) 的性质
DBL_MAX 在此代指浮点数能够表示的最大有限值。
我们之前提到了很多终止条件，他们都可能使得光线无法到达
且在落入
也就是说，如果规定 DBL_MAX，DBL_MAX，落入某一点时 DBL_MAX，那么无论如何
综上所述，我们可以通过
6. 二分法预处理
既然有 DBL_MAX， DBL_MAX，端点所属区间编号为
(1) 二分法使区间左右端点的区间编号与区间本身编号重合
对第一个区间进行二分法，得到区间中点
接着对第二个区间进行相同的处理，但此时不再更新第一个区间。这样处理完所有的区间后，所有的区间内都有一个解，且端点的区间编号都与其编号相同。
(2) 二分法使区间左右端点满足 Brent 法条件
接着对每个区间进行二分法，直到区间两端的 DBL_MAX 为止，此时已满足 Brent 法的条件，即可直接求解出第
7. 繁杂的数学推导
题解中我跳过了一些证明，如果你有兴趣，如下是详细证明。
(1) 从费马原理到光线方程
我们从一般的泛函开始
当其有
规定符号：
要注意，由于
由于存在约束，则对于变分
规定符号：
接下来，我们对辅助泛函进行变分，令函数
由于
由于起点终点位置固定，有
由此，我们得到了欧拉-拉格朗日方程与其边界条件，代入
这样我们就得到了有约束情况下的欧拉-拉格朗日方程。
从此出发，我们从费马原理求光线方程，对于费马原理：
我们可以知道
计算
那么对应的欧拉-拉格朗日方程为：
我们可以把中间的式子写作矢量式，通过给第
我们考虑沿着曲线
解上述方程，我们可以知道
于是我们得到了光线方程。
(2) 从光线方程到以 x 为自变量的微分方程
我们从光线方程出发，考虑空间维度
代入弧长
于是得到以
如果要以
解上述微分方程时，对于约束成立的初始条件，约束
(3) 从折射定律进行微元分析得到光子角动量守恒
不妨设光线随极坐标的
于是可知
实际上介质中光子动量为
在二维极坐标中，
最后附上代码
// 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;
}