边值问题是

$$\begin{cases} \dot{x}_i = \begin{cases} (0.5D^{-1}\psi)_i, \text{ if }(0.5D^{-1}\psi)_i \le 0 \\ 0 \text{, otherwise} \end{cases} \\ \dot{\psi} = 2\Sigma x \\ x(0) = x^0 \\ x(T) = 0 \end{cases}$$

这里$x^0>0$,$T$给出，$D>0$是对角矩阵，$\Sigma>0$是一个对称矩阵。实际上约束是$x(t)\ge 0\forall t$， 相当于$\dot{x}(t) \le 0$.

对于二维情况，我尝试用这段代码在 MATLAB 中解决它

function [tint, traj] = liq2mon_v3(x0, D, sigma, T)

  tlow = 0;
  thigh = T;
  L = length(x0); % L=2 in this example
  nsegm = 100; % I know it is for final calc only, not for the mesh
  solinit = bvpinit(linspace(tlow, thigh, nsegm), @guess);
  sol = bvp5c(@bvp4ode, @bvp4bc, solinit);
  tint = linspace(tlow, thigh);
  z = deval(sol, tint);
  traj = z(1:L, :);

  function xpsi = guess(t)
    xpsi = zeros(1, L*2);
    xpsi(1:L) = (x0*(1.0 - 1.0/T*t))';
  end

  function dydt = bvp4ode(t, y)
    xcomp = y(1:L); % L=2 in this example
    psicomp = y(L+1:end);
    u    = D \ psicomp / 2.0;
    dx   = min(u, 0); 
    dpsi = 2.0 * sigma * xcomp;
    dydt = [dx' dpsi']; 
  end

  function res = bvp4bc(ya, yb)
    xcompa = ya(1:2);
    xcompb = yb(1:2);
    res = [(xcompa - x0)' (xcompb)'];
  end

end

我得到一个错误

*Unable to solve the collocation equations -- a singular Jacobian  encountered
Error in ==> liq2mon_v3 at 13
sol = bvp5c(@bvp4ode, @bvp4bc, solinit)*

但是，如果我换行

dx = min(u, 0);

到

u(u>0) = 1e-13;
dx = u;

我确实得到了解决方案$\min(\min(x)) \approx -8\times 10^{-14}$对于相同的。$x^0$

问题：

1. 有没有办法在 MATLAB（或其他应用程序）中解决这个系统，以免违反约束？
   
2. 也许这可以在 Python 中解决，即使当前违反了约束？
3. 在这种情况下，为什么在 MATLAB 求解器中某处出现雅可比奇异数？