I am interested in computing a quantity that is mathematically defined as follows,

$$\phi(x_1,x_2) = \inf_{T>0} \inf_{\gamma \in C^{x_2}_{x_1 }(0,T)} \int_{0}^T L(\gamma, \dot{\gamma})\mathrm{d}t$$

where $x_1,x_2$ are points in $\mathbb{R}^n$ and $C^{x_2}_{x_1 }(0,T) = \{ \gamma: [0,T]\rightarrow \mathbb{R}^n| \gamma(0) = x_1, \gamma(T) = x_2 \}$ is the set of curves joining the two points $x_1, x_2$. $T$ represents, a time coordinate and the integrand is a an action $S$ with the corresponding Lagrangian. I am interested in paths which minimise the action, connecting $x_1, x_2$, and computing $\phi$ along those paths.

Supposing I know the Hamiltonian, one way I could think of doing this is writing the action in terms of: $S = \int_{0}^{t} \sum_i p_i\dot{q_i}-H(p,q)$, since essentially I am interested in curves in phase space that join two points and are parametrised by $T$, and minimise the action functional, they must satisfy Hamiltonian equations of motion (if I understand my mechanics correctly).

To solve the Hamiltonian equations of motion: $$\dot{q_i} = \partial_{p_i} H(p,q)$$ $$\dot{p_i} = -\partial_{q_i} H(p,q)$$

However, instead of initial conditions, I would like to specify a boundary value problem, $q(0) = q_1, q(T) = q_2, p(0) = p_1, p(0) = p_2$.

For this purpose I tried using NDSolve. Given a Hamiltonian, (j = 60, here $\rho,r$ stand in for $q, p$ respectively.

h2[j, \[Rho], 
  r] := ((-(1/900))*(-1 + 
      E^r)*(\[Rho]*(7921 + 7021*(-2 + \[Rho])*\[Rho]) + 
          E^r*(-1 + \[Rho])*(900 + 7021*\[Rho]^2)))/E^r

\[Rho]dt[j, \[Rho], 
  r] := ((1/900)*(\[Rho]*(-7921 - 7021*(-2 + \[Rho])*\[Rho]) - 
      E^(2*r)*(-1 + \[Rho])*(900 + 7021*\[Rho]^2)))/E^r

rdt[j, \[Rho], 
  r] := ((1/900)*(-1 + E^r)*(7921 + 7021*\[Rho]*(-4 + 3*\[Rho]) + 
      E^r*(900 + 7021*\[Rho]*(-2 + 3*\[Rho]))))/E^r

sp = StreamPlot[{rdt[j, \[Rho], r], \[Rho]dt[j, \[Rho], r]}, {r, -0.1,
    2}, {\[Rho], 0., 1}, StreamColorFunction -> "Purple", 
  PlotRange -> Full]

cp = ContourPlot[{rdt[j, \[Rho], r], \[Rho]dt[j, \[Rho], 
    r]}, {r, -0.1, 2}, {\[Rho], 0., 1}]
pp = ParametricPlot[{Log[(\[Rho] (-7921 - 
       7021 (-2 + \[Rho]) \[Rho]))/((-1 + \[Rho]) (900 + 
       7021 \[Rho]^2))], \[Rho]}, {\[Rho], 0, 1}, 
  PlotStyle -> {Thick, Red}]

enter image description here

Here, I have plotted the stream and contour plots for the two equations, along with the zero energy (H=0) contour. I would expect the minimum action plot between the the two fixed points, 0.15 and 0.5 ought to lie on the zero energy contour.

Now, I deploy NDSolve, following @Michael E2's second suggestion. (Choosing T to be sufficiently large, because techincally the transition is only possible in infinite T.

sols = FullSimplify[NDSolve[{D[[Rho][t], t] == D[h2[j, [Rho][t], r[t]], r[t]], D[r[t], t] == -D[h2[j, [Rho][t], r[t]], [Rho][t]], [Rho][0] == 0.151, [Rho][10^4] == 0.5}, {[Rho][t], r[t]}, {t, 0, 10^4}]]

Plotting the interpolating functions,

enter image description here

not quite what i wanted.

Outside of this particular example, is there built in functionality to find paths that minimize action in Mathematica, that can work for arbitrary systems in $\mathbb{R^n}$?

EDIT: It was recommended that I add some context to the problem. These types of optimisations appear in the context of stochastic dynamics and large deviation theory, and the quantity $\phi$ is the Friedlin-Wentzel quasipotential.

  • 1
    $\begingroup$ 1. For a system formed by two 1st order ODE, 2 constraints already determine a particular solution, if you give more constraints, the system will be overdermined. 2. Where's the definition of h2? $\endgroup$
    – xzczd
    Nov 12, 2019 at 4:41
  • 1
    $\begingroup$ You should give the definition of h2 in e.g. h2=… form instead of pasting the In[] and Out[], and avoid using TraditionalForm, these only make the code hard to check. Then, I'm not familiar with the topic discussed here, but is the following post related?: mathematica.stackexchange.com/q/191925/1871 $\endgroup$
    – xzczd
    Nov 12, 2019 at 5:09
  • 2
    $\begingroup$ The output can be transformed to InputForm, too. Just place the cursor at the output cell or select the output cell, and press Ctrl+Shift+I . $\endgroup$
    – xzczd
    Nov 12, 2019 at 6:53
  • 1
    $\begingroup$ Have you seen Variational Methods? I'm not sure if they'll help. -- I don't think solving the Hamiltonian equations of motion is the way to go. The integral curves in phase space do not go from a point $x_1$ to an arbitrary point $x_2$. (It is more or less equivalent to the NDSolve error you get.) But given a phase curve, it is the path of least action between the points it connects. $\endgroup$
    – Michael E2
    Nov 19, 2019 at 21:22
  • 1
    $\begingroup$ My approach, which I don't have time to explore, would be to have a number of variable points $p_j$, $j=1,\dots,n-1$, and fix $p_0=x_1$, $p_n=x_2$, and interpolate: Something line: gamma[p_?(MatrixQ[#, NumericQ]&), T_][t_] = Interpolation[Transpose[{Array[# &, {n+1, 1}, {{0., 1.}}], p}]][t/T], use NIntegrate to cmpute the action, and FindMinimum to minimize it. $\endgroup$
    – Michael E2
    Nov 20, 2019 at 13:35


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.