Dynamical Systems in Lean4

Stream: Natural sciences

[Himanshu Jain]

Hi everyone,

I have currently enrolled in a Lean 4 Theorem Prover course in my institute. As part of the course, I need to choose a project.
I’m considering formalizing the derivation of Hamilton's equation of motion from the Lagrangian formulation (in the way presented in Goldstein's Classical Mechanics Chapter 8 The Hamilton Equations of Motion). I am not sure whether we can do this to full generality because then this might become the problem in symbolic computation also. Maybe we have to take some system say Harmonic oscillator and then from the Lagrangian, we have to construct the Hamiltonian (Legendre transform of Lagrangian) and then proceed further.

From my search, I found the scilib library, which includes theorems on derivatives and calculus, but it seems to lack results on multivariable calculus.
Given that this is a 1-month project, I’d like to know how feasible this would be. Are there any major challenges I should expect?
Also, if you have any suggestions for alternative project topics, I’d love to hear them!

Thanks in advance.

[Tomas Skrivan]

I think this should be possible and I'm willing to help. To get you started I wrote down the formal statement of Lagrange and Hamilton's equations:

import Mathlib

variable {X} [NormedAddCommGroup X] [InnerProductSpace ℝ X] [CompleteSpace X]

open RealInnerProductSpace

def LagrangeEquations (L : X → X → ℝ) (x : ℝ → X) : Prop :=
  let v := deriv x
  ∀ t, deriv (fun t' => gradient (fun v' => L (x t') v') (v t')) t -- d/dt (∂L/∂ẋ)
       -
       gradient (fun x' => L x' (v t)) (x t)                      -- ∂L/∂x
       =
       0

def HamiltonEquations (H : X → X → ℝ) (q p : ℝ → X) : Prop :=
  ∀ t, deriv q t =   gradient (fun p' => H (q t) p') (p t)
       ∧
       deriv p t = - gradient (fun x' => H x' (p t)) (q t)

As you can see, there is no need for explicit partial derivative/multivariate calculus. You just work with Lagrangian and Hamiltonian as a function of two argument X → X → ℝ and take corresponding derivatives w.r.t. the first or the second variable which effectively gives you partial derivatives.

The function that obtains Hamiltonian from Lagrangian is:

noncomputable
def LtoH (L : X → X → ℝ) (q p : X) : ℝ :=
  let v := (gradient (fun v' => L q v')).invFun p
  ⟪p, v⟫ - L q v

where let v := (gradient (fun v' => L q v')).invFun p is inversion of the expression $p = \frac{\partial L}{\partial \dot q}$ which is at the heart of Lagrange transform.

Now we can state the main result, if a trajectory x : ℝ → X satisfies Lagrange equation then it satisfies Hamilton's equations

theorem LE_implies_HE (L : X → X → ℝ) (x : ℝ → X)
  (h : LagrangeEquations L x) :
  let q := x
  let p := fun t => gradient (fun v' => L (x t) v') (deriv x t)
  HamiltonEquations (LtoH L) q p := sorry

your task is to fill in the sorry. It is not provable as stated, you will need some technical conditions on L and x such as appropriate differentiability and some form of convexity of L in v argument.

It would be nice to have the result also in revese:

noncomputable
def HtoL (H : X → X → ℝ) (x v : X) : ℝ := sorry

theorem HE_implies_LE (H : X → X → ℝ) (p q : ℝ → X)
  (h : LagrangeEquations L x) :
  LagrangeEquations (HtoL H) q := sorry

where HtoL is just analogous Legendre transform taking you from Hamiltonian to Lagrangian.

---

I'm expecting that the hardest part will be reasoning about the function gradient (fun v' => L q v'). Under what conditions on L is it invertible? Under what conditions is the inverse differentiable and what is its derivative?

I think a reasonable strategy would be to just assume all the steps involving inverse function, like (f⁻¹)'(y) = 1/f'(f⁻¹(y)), and get the main structure of the proof. Once that is done, you can maybe prove everything formally correct in case of X = ℝ. Doing it for general space X might be a bit challenging as, I think, the necessary theorems are not formulated in the most beginner friendly way.

[Tomas Skrivan]

This would also make a nice PR to PhysLean, right? @Joseph Tooby-Smith

[Tomas Skrivan]

Also, no matter what I would start with X = ℝ. In this case, gradient can be replaced with deriv. The reason is that gradient is relatively new and I'm not sure if the API is well developed thus might be hard for a beginner to use.

[Joseph Tooby-Smith]

Tomas Skrivan said:

This would also make a nice PR to PhysLean, right? Joseph Tooby-Smith

Indeed this would be a very nice PR to PhysLean. I actually started some basic stuff in the direction of classical mechanics yesterday, and indeed I did stuff related harmonic oscillator. See here. I only did the basics though: I wrote down the lagrangian and the solutions but did not connect the two with the Euler-Larange equations. So anything in this direction would be great :). I don't really plan to work more in the direction of classical mechanics for the time being.

[Himanshu Jain]

@Tomas Skrivan , @Joseph Tooby-Smith , Thanks a lot for your detailed response!

I come from a scientific computing background, where I primarily work with numerical methods like Finite Differences to solve PDEs. When I first started learning Lean4 at my institute, I struggled with its proof-oriented paradigm—perhaps because I’ve been working with Python for so long.
When I first thought about formalizing physical concepts in Lean for my course project, I found it difficult to visualize how this would work. Initially, I assumed we might define a specific system, express the Lagrangian and Lagrange's equation symbolically (like in SymPy), and manipulate expressions to derive Hamilton’s equations. However, I wasn’t even clear on what the actual inputs and outputs of our Lean functions would be.

From this discussion, I now understand that:

• Input: Formal definitions of Lagrange’s equations, Hamilton’s equations, and the Legendre transform.
• Output: A proof that, under certain conditions, if a function $x(t)$ satisfies Lagrange's equations, then the corresponding $q(t)$ and $p(t)$ satisfy Hamilton’s equations.

Of course, we must assume constraints on the interaction potential of the system to ensure the Lagrangian is invertible. Does this interpretation align with what you had in mind?

The project seems quite interesting, and I’d love to discuss it further with you to gain better insights!

[Joseph Tooby-Smith]

This does indeed seem like an interesting project :).

I agree with your input and output, although I would be careful with the phrase Lean functions, as it is not like Lean has an explicit function that takes these inputs and gives this output.

Also, in case you haven't seen it there is Zulip channel for PhysLean here.

[Tomas Skrivan]

I'm not sure if I would call it input and output. In a way everything is an input, you have to write the statement and the proof, and the only thing Lean outputs is yay or nay based on the correctness of your proof.

[Alex Meiburg]

Well. There could be something to be said for actually writing "implementations" too. :) I mean, Lean does generate code. But that would be exceedingly hard to do with derivatives, real numbers, etc. all being noncomputable.

[Alex Meiburg]

Tomas Skrivan said:

I'm expecting that the hardest part will be reasoning about the function gradient (fun v' => L q v'). Under what conditions on L is it invertible? Under what conditions is the inverse differentiable and what is its derivative?

I think this part is worth highlighting extra. A "typical" treatment of Hamiltonian/Lagrangian mechanics at the level of Goldstein's sort of throws out many of these assumptions. You're really working with the Legendre transform here, which will require some assumptions about convexity and invertibility and so on if you want to get the proofs to check.

I think that, before drafting the proofs in Lean, you should look at the book proofs you want to formalize, and think quite critically about what assumptions are needed on these solutions at each step.

[Alex Meiburg]

For instance: there's a classic example where a ball starts at rest atop a cone shaped like y = -x^(4/3). According to Hamilton's equations, there's an infinite number of solutions - at any time t, it could 'spontaneously' start rolling down the cone. (And in reverse, it could have been rolling up the cone and come to rest at any point in time in the past.) So, despite being a pretty nice differentiable potential, there is no uniqueness on the differential equations. You should think about how you'll handle that!

[Tomas Skrivan]

Alex Meiburg said:

I think that, before drafting the proofs in Lean, you should look at the book proofs you want to formalize, and think quite critically about what assumptions are needed on these solutions at each step.

I think the exact opposite of this approach should be doable too. I would write the proof the way you find it in standard physics textbook first and sorry all the assumptions of used theorems. After that I would inspect all the sorry, realized what assumptions do I need on the potential to allow me to remove every sorry.

[Tomas Skrivan]

Alex Meiburg said:

For instance: there's a classic example where a ball starts at rest atop a cone shaped like y = -x^(4/3). According to Hamilton's equations, there's an infinite number of solutions - at any time t, it could 'spontaneously' start rolling down the cone. (And in reverse, it could have been rolling up the cone and come to rest at any point in time in the past.) So, despite being a pretty nice differentiable potential, there is no uniqueness on the differential equations. You should think about how you'll handle that!

Yes this is indeed an issue but it is not related to the original problem:

I’m considering formalizing the derivation of Hamilton's equation of motion from the Lagrangian formulation

Which, in my interpretation, boils down to showing that if a trajectory satisfies Lagrange equation then it satisfies Hamilton's equations. This is unrelated to the question about existence and uniqueness as you are given a solution as a hypothesis.

I'm expecting that the only requirement will be:

• smoothness of Lagrangian ContDiff ℝ 2 (fun (x,v) => L x v))
• strong convexity, expressed as positive definitness of the hessian:

∃ c > (0:ℝ), ∀ (x v dv : X), fderiv ℝ (fderiv ℝ (L x ·) · dv) v dv ≥ c * ‖dv‖^2

I'm not sure if mathlib has theorems about functions that satisfy this kind of convexity and about positive definite bilinear forms. That is why I'm suggesting to omit those proofs and maybe prove them only in a special case X = ℝ where the convexity condition just says that iteratedDeriv 2 (L x ·) v > 0.

---

If we want to talk about local existence and uniqueness of the solution we can use Picard–Lindelöf theorem which requires that the equation is Lipschitz in space. In your example the energy is $- |x|^{\frac43}`$ thus the force is $- \frac43 \frac{x}{|x|^\frac23}$ which is not Liptschitz thus not satisfying Picard–Lindelöf theorem.

Under the ContDiff ℝ 2 (fun (x,v) => L x v)) assumption you get Liptschitz condition immediately thus you are guaranteed local existence and uniqueness.

Global existence would also require some integrability condition which I think is not yet in mathlib.

[Tomas Skrivan]

Alex Meiburg said:

Well. There could be something to be said for actually writing "implementations" too. :) I mean, Lean does generate code. But that would be exceedingly hard to do with derivatives, real numbers, etc. all being noncomputable.

This is what my project, SciLean, is all about. Have a look at here for generating a numerical program from just specifying Hamiltonian.