I am trying to implement an initial problem. Let the system of ODEs be

\frac{d\vec{y}}{dt} = A(t;\vec{\theta})\vec{y}, \quad \vec{y}(0) = \vec{y}_0, \quad 0 \leq t \leq T.

The matrix A depends on the time t, as well a vector of control parameters \vec{\theta}.

The time-stepping method I’m using is similar to the trapezoidal rule. The timestepping is done by something like

(I - A(t_{n+1}, \vec\theta))\vec y_{n+1} = (I + A(t_n))\vec y_n

For reasons particular to my use-case, I would like the linear solves in my timestepping method to be iteratively using `gmres!`

from the `IterativeSolvers.jl`

package. What’s more, I would like to compute d\vec{y}/dt without actually computing the matrix A. I have a function `apply_A!(dydt, y, t, theta)`

, which compute the derivative of `y`

and stores it in `dydt`

.

I have been trying to do this using `LinearMaps.jl`

, in a way like this:

```
using LinearMaps
using IterativeSolvers
#=
The problem being solved is
y' = -2t*y
y(t) = y₀*exp(-t²)
The scalar equation is done for each element of y, so you can solve the initial
value problem for multiple initial conditions.
=#
function apply_A!(dydt, y_in, t)
dydt .= (-2*t) .* y_in
return nothing
end
function solver(y0, dt, nsteps)
y = copy(y0)
dydt = zeros(length(y))
t = 0.0
right_hand_side = zeros(length(y))
function compute_left_hand_side!(out_vec, y_in)
apply_A!(dydt, y_in, t)
out_vec .= y_in
out_vec .-= 0.5*dt .* dydt
return nothing
end
LHS_map = LinearMaps.LinearMap(
compute_left_hand_side!,
length(y), length(y),
ismutating=true
)
for n in 0:(nsteps-1)
t = n*dt
apply_A!(dydt, y, t)
right_hand_side .= y
right_hand_side .+= 0.5*dt .* dydt
#display(y)
t = (n+1)*dt
gmres!(y, LHS_map, right_hand_side)
end
println("Final t = ", t)
println("Final y = ", y)
return y
end
```

I have something like this implemented, but my method is pretty slow, makes many more memory allocations than are reasonable, and I have tracked down that the vast majority of the memory allocations do not come from any of the code I have written (i.e. they are from gmres or the linear map, in the way that I have used them).

I am not very experienced with debugging type-allocation, but it seemed to me from the beginning that this might be a bad idea, since I create a function inside a function, which I then pass to `gmres!`

. And so the function `gmres!`

takes may not be known until runtime, messing with the compilation process. Also, `compute_left_hand_side`

relies on variables outside the scope of that function. Although those variables are technically not global variables, which I know people generally advise against using, I fear it may be similarly bad for performance to use variables outside the function scope like I have here.

Does anyone have a “right” way to do this? I feel like I am working against the intended use of these packages, but I can’t think of another way to accomplish this. If there are packages besides `LinearMaps`

and `IterativeSolvers`

which would be more appropriate for this use-case, I am happy to switch.

# EDIT

Changed code mock-up to a MWE (and one that is actually algorithmically correct).