Taylor Series Approximations

A common issue when trying to verify interesting dynamics in KeYmaera X is that dL and KeYmaera X do not provide any native support for trigonometric or other transcendental functions, and thus such functions must be expressed as differential equations instead. There is a good reason for this: arithmetic problems over polynomials are decidable, but over transcendental functions they are not.

If we wish to reason about transcendental functions using the polynomials that KeYmaera X knows and loves, we need to approximate them from above or below. Luckily, there is a general approach for upper or lower-bounding arbitrary (continuous) functions: Taylor Series! Recall that a function can be rewritten as a power series using its derivatives:

If we take the first k terms of this series for some constant k, we will not recover f(x) exactly, but it is generally easy to tell whether the first k terms are an upper bound vs. a lower bound. In either case, we can prove that this bound is in fact a bound, which can often be used in place of an exact solution when implementing a controller + proving a safety theorem.

Example 1: Natural logarithm ln(1 - x)

The function f(x) = ln(1 - x) has the following Taylor series when 0 < x < 1:

And since ln(1-x) = 1/(x-1), we express it with the ODE {x' = 1, ln' = 1/(x-1)}. The following KeYmaera X model expresses the 3rd-degree taylor bound for ln(1-x):

ProgramVariables
  Real x;
  Real ln;
End.

Problem
  x = 0 & ln = 0 -> [{x' = 1, ln' = 1/(x-1) & x>=0&x<1}] ln <= -x - x^2/2 - x^3/3
End.

The strategy for proving Taylor bounds is simple: progressively cut in all your lower-degree Taylor bounds, proving each bound (including the final one) by dI. Recall that as a shortcut we can use diffInvariant to cut and prove an invariant all in one step. Thus the following tactic will prove the above bound:

unfold;
diffInvariant("2/(x-1)^3<=-2", 1);
diffInvariant("-1/(x-1)^2<=-1-2*x", 1);
diffInvariant("1/(x-1)<=-(1)-x-x^2", 1);
dI(1)

Example 2: Hyperbolic cosine and sin cosh and sinh

Recall the functions cosh and sinh (hyperbolic cosine and sin). They are nice and easy to express as an ODE, because they satisfy the following equations:

cosh'(x) = sinh(x)
sinh'(x) = cosh(x)
cosh(0) = 1
sinh(0) = 0

Thus the following ODE solves to cosh = cosh(t), sinh = sinh(t):

cosh = 1 & sinh = 0 & x = 0 -> [{cosh' = sinh, sinh' = cosh, x' = 1}]...

How do we know that this ODE solves exactly to cosh(x) and sinh(x) instead of some other pair of functions? Recall that ODE solutions (under the right conditions, which are always satisfied in DL) are unique. The identities above say that cosh and sinh must be a solution to the ODE, and thus they are the only solution.

cosh and sinh also have nice Taylor series representations:

Which is to say that cosh(x) + sinh(x) = e^x, and cosh gets all the even terms of the e^x power series while sinh gets the odd ones. The following KeYmaeara X model expresses the third Taylor bound for sinh:

ProgramVariables
  Real x;
  Real cosh;
  Real sinh;
End.

Problem
  x = 0 & cosh = 1 & sinh = 0 ->
  [{cosh' = sinh, sinh' = cosh, x' = 1 & cosh >= 1}]
  sinh >= x + (x^3)/6 + (x^5)/120
End.

This ODE is a bit more complicated than that for ln, because it has two mutually-recursive variables. You might wonder if we'd get stuck in a loop trying to prove a Taylor bound for a recursive ODE. Thankfully we don't because the order of our Taylor bound decreases at each step. The idea is that to prove the order-k bound for cosh we cut in the order-(k-1) bound for sinh and vice-versa. Thus the following series of cuts will prove the above model:

unfold;
diffInvariant("sinh >= x", 1);
diffInvariant("cosh >= 1 + (x^2)/2", 1);
diffInvariant("sinh >= x + (x^3)/6", 1);
diffInvariant("cosh >= 1 + (x^2)/2 + (x^4)/24", 1);
dI(1)

Note that the first cut relies on the domain constraint cosh >= 1 which I have conveniently added to my model. This serves as the "base case" for the Taylor bound process. Formula cosh >= 1 can also be added as cut, but this might require you to be a lot more clever since you will probably need an invariant that talks about both cosh and sinh to prove sinh >= 1. Therefore when working on a proof it's best to start out with making cosh >= 1 an assumption in your domain constraint and then once you've verified the model, adapt it so you prove cosh >= 1 as a cut instead of making it a constraint in the model.