计算科学 - 在 Matlab 中解决一个明显棘手的测地线 BVP - 吾爱随笔录

我希望能够解决 BVP

{\ddot{μ}}_{k} = - \frac{μ_{k}}{2} (\sum_{i = 1}^{n} \frac{{\dot{μ}}_{i}^{2}}{μ_{i}} - \frac{{\dot{μ}}_{k}^{2}}{μ_{k}^{2}} + \frac{{[\sum_{i = 1}^{n} {\dot{μ}}_{i}]}^{2}}{μ_{0}})

$\ddot \mu_k = -\frac{\mu_k}{2} \left ( \sum_{i=1}^n \frac{\dot \mu_i^2}{\mu_i} - \frac{\dot \mu_k^2}{\mu_k^2} + \frac{ \left [ \sum_{i=1}^n \dot \mu_i \right ]^2}{\mu_0} \right )$ 和

μ_{0} := 1 - \sum_{i = 1}^{n} μ_{i}

$\mu_0 := 1-\sum_{i=1}^n \mu_i$ ,

μ_{j} \in [0, 1]

$\mu_j \in [0,1]$ 同样对于

0 \leq j \leq n

$0 \le j \le n$ ，并且具有受上述约束的任意边界值。因为这个 BVP 是某个黎曼（伪）流形的测地线方程，所以我知道它每次都有唯一的解。事实上，我知道解析的解决方案是什么（基本上，坐标

\sqrt{μ_{i}}

$\sqrt{\mu_i}$ 跟随大圆圈），但我真正关心的是获得与此一致的数值结果，以便我可以概括问题并仍然获得数值（没有分析结果的希望）。

但是，我一生都无法弄清楚如何让搭配方程在 MATLAB 中工作，我怀疑在这种情况下拍摄方法是否明智 $n>1$ （无论如何我都不知道如何让它们在这种情况下工作，尽管我已经得到了一个 IVP 似乎不错的结果 $n = 1$ ）。

这是我到目前为止所得到的：

%% Set up symbolic system of ODEs
n = 1;
eval(['syms mu(t) [1,',num2str(n),']']);
mu0 = 1-sum(mu); 
geodesicEqn = diff(mu,2) == ...
    -(mu/2).*(sum((diff(mu,1).^2)./mu)...
    -((diff(mu,1).^2)./(mu.^2))...
    +(sum(diff(mu,1))^2)./mu0);   % (sum(diff(mu,1))^2)./mu0 = (diff(mu0,1).^2)./mu0

%% Translate ODE system to vector field and MATLAB function in turn
V = odeToVectorField(geodesicEqn)
M = matlabFunction(V,'vars',{'t','Y'})

%% Set up random but reproducible initial conditions
interval = [0,1];
rng('default');
init = rand(1,numel(V)/2)/(numel(V)/2)
term = rand(1,numel(V)/2)/(numel(V)/2)
a = [init,1-sum(init)];
b = ([term,1-sum(term)]-a)*interval(end)+a;

%% Boundary value problem (FAILS FOR EVERY GUESS I'VE TRIED)
guess = @(t) [cos(pi*t);sin(pi*t)];
solinit = bvpinit(linspace(eps^.125,1-eps^.125,10),guess);
bcfun = @(y0,y1) [y0(1)-init;y1(1)-term];
bvpopts = bvpset('NMax',1e5);
sol = bvp4c(M,bcfun,solinit,bvpopts)

无论我做出什么猜测，MATLAB 都会抱怨搭配方程。任何提示都会很棒。