High precision square-root integration

[andreasvarga]

This problem arises in the larger context of solving differential matrix Lyapunov equations with periodic coefficients, where the solution is guaranteed to be nonnegative definite. Therefore, it would be desirable to solve these equations directly for their Cholesky factors, to gain additional accuracy which is typically provided by working with square-root quantities. The following example is (probably) the simplest case which can be imagined (first order, constant coefficients):

dx/dt = -x +1, x(0) = 0, 0 <= t <= 1

which has x(t) = 1-exp(-t) as solution.

If we express x = u^2, the square-root quantity u satisfies the implicit differential equation

2*u*du/dt = -u^2 + 1, u(0) = 0, 0 <= t <= 1

The following code fails to compute the solution

using DifferentialEquations
using Sundials
function f(out,du,u,p,t)
  out[1] = -2*u[1]*du[1]-u[1]^2+1
end
tspan = (0.0, 1.0)
differential_vars = [true]
u0 = [0.]
du0 = [0]
prob = DAEProblem(f, du0, u0, tspan, differential_vars = differential_vars)
sol = solve(prob, IDA(), abstol=1.e-10,reltol=1.e-10)

with the error

[IDAS ERROR]  IDACalcIC
  Newton/Linesearch algorithm failed to converge.

retcode: InitialFailure
Interpolation: 3rd order Hermite
t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{Vector{Float64}}:
 [0.0]

Perturbing the initial condition to
u0 = [0.00001]
a solution can be computed, which is however far from a "high precision" solution

julia> sqrt.(1 .-exp.(-sol.t)) - vcat(sol.(sol.t)...)
1141-element Vector{Float64}:
 -1.0e-5
 -9.96843986386683e-6
 -9.955396516668458e-6
 -9.93697972687653e-6
 -9.910993028213746e-6
 -9.874359429691951e-6
  ⋮
  7.506912869104099e-10
  7.810937452390476e-10
  7.887508424175849e-10
  8.013723018507335e-10
  8.214686708640784e-10

I wonder if there is a way to solve the implicit differential equation to high accuracy on the whole time span. Thanks in advance for any hints.

[andreasvarga]

For this simple example I can determine a high precision solution using a consistent initialization of the derivative. With the following code

using DifferentialEquations
using Sundials
function f(out,du,u,p,t)
  out[1] = -2*u[1]*du[1]-u[1]^2+1
end
tspan = (0.0, 1.0)
differential_vars = [true]
u0 = [1.e-11]
du0 = [1/u0[1]/2]
prob = DAEProblem(f, du0, u0, tspan, differential_vars = differential_vars)
sol = solve(prob, IDA(), abstol=1.e-10,reltol=1.e-10)

I got the following results

julia> sqrt.(1 .-exp.(-sol.t)) - vcat(sol.(sol.t)...)
1326-element Vector{Float64}:
 -1.0e-11
 -4.175143073132806e-11
 -6.47703804423539e-11
 -8.817497863579657e-11
 -8.938104423951714e-11
 -9.134777588272761e-11
  ⋮
  7.055517281528978e-10
  7.05854708016318e-10
  7.058055251363271e-10
  7.054212769475043e-10
  7.019628212034945e-10

which seems to be satisfactory.

Just a final remark for those interested in a possibly open research issue:

The next step for me would be the extension to the matrix case x = u'*u, with u upper triangular, for which
dx/dt = (du/dt)'*u+u'*(du/dt)
and a consistent initialization would mean to solve a Lyapunov-like matrix equation satisfied by du0
du0'*u0+u0'*du0 = q
where q is a symmetric nonnegative definite matrix and du0 is upper triangular. See also the discussion at #23.