DIfferentialEquations with the driving function defined on a fixed grid

General Usage
1

Hello,

I want to use DifferentialEquations.jl to solve a differential equation whose driving function is defined on a fixed grid. One of thе solutions would be to pass the grid and the driving function arrays as parameters and then do an interpolation at every integration step. Something like this:

using DifferentialEquations
using PyPlot

alg = BS3()


# analytic driving function ----------------------------------------------------
function func1(u, p, t)
    return cos(t)
end

tspan = (0.0, 10.0)
u0 = 0.0

prob1 = ODEProblem(func1, u0, tspan)
sol1 = solve(prob1, alg, saveat=0.1)


# driving function defined on the grid -----------------------------------------
function func2(u, p, t)
    tt, ff = p

    # linear interpolation:
    if t <= tt[1]
        f = ff[1]
    elseif t >= tt[end]
        f = ff[end]   
    else
        ndt = (t - tt[1]) / (tt[2] - tt[1])   # number of steps from tt[1] to t
        i = Int(floor(ndt)) + 1
        f = ff[i] + (ff[i+1] - ff[i]) / (tt[i+1] - tt[i]) * (t - tt[i])
    end

    return f
end

t = range(0.0, 10.0, length=101)
f = cos.(t)

tspan = (t[1], t[end])
u0 = 0.
p = (t, f)

prob2 = ODEProblem(func2, u0, tspan, p)
sol2 = solve(prob2, alg, saveat=(t[2]-t[1]))


# ------------------------------------------------------------------------------
plot(t, f)   # driving function
plot(t, sin.(t))   # exact solution
plot(sol1.t, sol1.u)
plot(sol2.t, sol2.u)
show()

However, are there any more smart ways to solve such problems?

2

Certainly you will have to do some kind of interpolation (linear, cubic, etc.), e.g. via Interpolations.jl or Dierckx.jl. However, beware that such interpolations have discontinuities in some derivatives at the “knots” of the interpolant (typically knots = the data points or a subset thereof). This will greatly degrade the accuracy/efficiency of the solver algorithms in DifferentialEquations.jl, since they assume smoothness. A workaround is to pass a tstops parameter listing all of the knots, which will allow the solver to integrate over the smooth segments.

If you have a lot of knots (finely spaced data), however, breaking integration at every knot might cause the solver to use more steps than it would for a similar ODE with smooth data.

An alternative is to find an interpolation that is smooth everywhere, but this is hard to do in general. Simply fitting to a high-degree polynomial is likely to suffer from a Runge phenomenon in the common case of equally spaced data points. If your data represents a smooth function and if you can choose the grid of data points, then by choosing a set of Chebyshev points you can use a smooth Chebyshev interpolant; ApproxFun.jl may help with this. Alternatively, if your data are noisy, then you might consider some kind of smooth fit.

(In the special case where you are just computing a definite integral, i.e. a separable ODE, then you can use QuadGK.jl to design a specialized quadrature rule that takes your nonsmooth interpolant into account a priori.)

5 Likes
3

Thank you for the detailed answer. I did not expect that there are so much stuff underneath it all. Can you recommend something more to read about it?

Meanwhile, assuming that the grid is fine enough, what is the way to force the solver to step over only the predefined steps? I want to be sure, that under the hood, the solver does not split the step, independently of the chosen algorithm. The documentation mentions the fixed step size only for the methods without adaptivity (https://docs.juliadiffeq.org/latest/basics/common_solver_opts/#Fixed-Stepsize-Usage-1).