How to set up the differential equation system to speed up computation?

Question

I've set up a system of differential equations, obtained after discretizing pde, in the following way

dY(1) =  terms_starty;
dY(end) = terms_endy;
dY(2:end-1) = (1./Vector1).*(Matrix1*Y+Matrix2*Y);

dZ(1) =  terms_startz;
dZ(end) = terms_endz;
dZ(2:end-1) = (1./Vector1).*(Matrix1*Z+Matrix2*Z);

dYZ = [dY dZ];

I'm solving the above system in MATLAB using ode15s.

odeset('abstol', 1e-10, 'reltol', 1e-9)
[t, s]  = ode15s(@(t,s) factory(t,s), tspan , s0, options);

The vector Y,Z is of size ~1200 and the number of differential equations is ~2400

The total simulation time takes around 570s when

odeset('abstol', 1e-10, 'reltol', 1e-9)

is used for error settings and it reduces by 5 times when the default error settings are used. I tried to use the profiler and the ode solver takes around 508 s.

I'd like to know if there are ways using which I can speed up the compute time taken by the ode solver. Also, could someone kindly clarify if the ode solver calls BLAS functions for Matrix operations? I'm using version 2019b of MATLAB.

EDIT: I'm trying to understand how to generate the jpattern matrix with the simple example below

  syms x y z;
    F = [x*y, cos(x*z), log(3*x*z*y)]
    v = [x y z]
    J = jacobian(F,v)

gives,

  J =
     
    [           y,   x,           0]
    [ -z*sin(x*z),   0, -x*sin(x*z)]
    [         1/x, 1/y,         1/z]

From J I'd like to generate the jpattern matrix of the following form,

jpattern = 
[ 1,   1,  0]
[ 1,   0,  1]
[ 1,   1,  1]

But I couldn't find the appropriate command that I could use to convert J to the form in jpattern. Basically, I'm trying to find out a way to assign 1 in place of the non-zero entries in J. Suggestions on how to do this will be really helpful.

EDIT2: I tried to set up the sparsity pattern in the toy model below,

x0 = [1 0 0 0 0 0 0 0 0 0]';
tspan = 0:0.01:5;

f0 = fun(0, x0);
J = odenumjac(fun,{0 x0}, f0); 
sparsity_pattern = sparse(J~=0.);
options = odeset('Stats', 'on', 'JPattern', sparsity_pattern);
[t, sol]  =  ode15s(@(t,x) fun(t,x), tspan , x0);
plot(t, sol)

function f = fun(t,x)

mat1=[ 
       1    -2     1     0     0     0     0     0     0     0;
       0     1    -2     1     0     0     0     0     0     0;
       0     0     1    -2     1     0     0     0     0     0;
       0     0     0     1    -2     1     0     0     0     0;
       0     0     0     0     1    -2     1     0     0     0;
       0     0     0     0     0     1    -2     1     0     0;
       0     0     0     0     0     0     1    -2     1     0;
       0     0     0     0     0     0     0     1    -2     1;
       ];

mat2 = [
        1    -1     0     0     0     0     0     0     0     0;
        0     1    -1     0     0     0     0     0     0     0;
        0     0     1    -1     0     0     0     0     0     0;
        0     0     0     1    -1     0     0     0     0     0;
        0     0     0     0     1    -1     0     0     0     0;
        0     0     0     0     0     1    -1     0     0     0;
        0     0     0     0     0     0     1    -1     0     0;
        0     0     0     0     0     0     0     1    -1     0;
        ];

f(1,1) = 0;     
f(2:9,1) = mat1*x + mat2*x;
f(10,1) = 2*(x(end-1) - x(end));
end

Unfortunately, I couldn't figure out how to fix the error

Not enough input arguments.

Error in Untitled>fun (line 36)
f(2:9,1) = mat1*x + mat2*x;

Error in Untitled (line 5)
J = odenumjac(fun,{0 x0}, f0);

Suggestions on how to fix this error will be of great help.

I wasn't sure if the vectorization was done right. So I tried the following as well

function df = fun(t,x)
  global mat1 mat2
  f(1,:) = 0;
  f(2:9,:) = mat1*x + mat2*x;
  f(10,:) = 2*(x(end-1) - x(end));
  df = [f(1, :); f(2:9, :); f(10, :)];
end

This didn't help me either.

Answer 1

Some things I can think of:

• use sparse matrices for Matrix1 and Matrix2 to speed up the computations of dZand dY
• use larger integration tolerances reltoland abstol, especially if your are searching for steady-state solutions and/or do not need a precise resolution of the transient dynamics of your system.
• you are solving with ode15s, which is an implicit integrator, you should provide the Jacobian of your ODE function yourself, as it seems pretty straightforward to do here. This is done with the keyword "Jacobian".
• if the Jacobian is too heavy to derive (maybe for a more complex problem), you can at least specify its sparsity pattern via the keyword "JPattern". If the jacobian is block-diagonal, the solver will be able to compute it via finite-differences in much fewer calls than it would considering a dense Jacobian (which is the default otherwise).
• finally, vectorise your function and set the keyword "Vectorized" to True. This way, the solver can call your ODE functions with multiple inputs at a time and sparse some overhead. This is especially useful for the previous Jacobian computation.

Here, I think the code can easily be vectorized by simply rewriting your equalities to (+ or - some adjustments for the element-wise operations):

dY(1,:) = ...
dY(end,:) = ...
dY(2:end-1,:) = ... 

Oh and of course, if your problem is a non-stiff ODE, just use explicit solvers like ode45, they will be much quicker. Unless you want to achieve steady-state and the initial transient is a long, in which case implicit method with large error tolerances may be able to reach steady-state much more quickly, without having to finely resolve the transient.

EDIT: if I misunderstood your answer and you actually want to run 2400 different simulations, it may be a bad idea to run them as a single "big" smulation, as the time step will be limited by the "hardest" solution. An alternative, in that case, is to use a parallel for loop (using parfor) and to solve each case separately over multiple cores.

EDIT 15/02/2021 Regarding odenumjac and the formation fo the sparsity apttern, here is how you can do it:

f0 = ode(0,y0); % value of your ODE function
J = odenumjac(ode,{0 y0}, f0); ! Computes the Jacobian, assuming it is dense
sparsity_pattern = sparse(Jac~=0.);

You can then specify 'JPattern', sparsity_pattern as an additional option when calling the ODE solver. Matlab will take advantage of the knowledge of this pattern to approximate the Jacobian faster via finite-differences. Note that, in your case, your Jacobian is easy to derive, so you would be better off providing your Jacobian through a function via the option Jacobian or, if it is constant, as a constant Jacobian via JConstant.

Regarding the way odenumjac approxiamtes the Jacobian, your best bet is to lookup "Jacobian finite-difference" on the internet to find more in-depth explanations, but here is the main principle. You ODE function $f$ takes an input (state vector $x$) of size $N$. Some ODE solvers (usually the implicit ones), require the knowledge of the Jacobian $\partial_x f` ($NxN$ matrix) to be able to compute a solution.

The $i$-th row of $\partial_x f$ is the vector of the partial derivatives of the $i$-th component of $f$ with respect to each component of $x$. Conversely, $i$-th column of $\partial_x f$ is the vector of the partial derivatives of all the components of $f$ with respect to the $i$-th component of $y$. This column can be approximated via the following first-order finite difference formula:

$$(\partial_{x_i} f)(x) = \dfrac{f(x+ h e_i) - f(x)}{h}$$ with $e_i = (0,..,0,1,0,..0)$,the vector with all components set to zero excapt the $i$-th one, which is set to one.

Then, doing a loop on $i$, the whole Jacobian can be approximated. This would required $N+1$ evaluations of $f$: one unperturbed to compute $f(x)$, and $N$ evaluations to compute $f(x+h e_i), i=1..N$.

For general problems with a sparse Jacobian structure (for example discretised PDEs), the knowledge of the sparsity pattern can be used to perform multiple perturbations at once, when different components of $x$ do not have an "overlapping" impact on the components of $f$. You can see the article "On the estimation of sparse Jacobian matrices" for more information: http://www.ii.uib.no/forskningsgrupper/opt/forskning/papers/extern/cpr.pdf

Answer 2

The simple example as posted wasn't stiff, but I put it together in Julia anyways for show. I modified it to be a system of 2 PDEs with N=1200 and get in the 10's of ms

using ModelingToolkit, LinearAlgebra, BenchmarkTools, DifferentialEquations

# Setup matrices
N = 1200
mat1 = hcat(zeros(N),Tridiagonal(ones(N-1),-2*ones(N),ones(N-1)),zeros(N))
mat1[1,1] = 1
mat1[end,end] = 1

mat2 = hcat(zeros(N),Bidiagonal(-1*ones(N),ones(N-1), :L),zeros(N))
mat2[1,1] = 1

# Define the ODE

function f(du,u,p,t)
    y = @view u[:,1]
    z = @view u[:,2]
    dy = @view du[:,1]
    dz = @view du[:,2]

    dy[1] = 0
    dy[2:end-1] .= mat1*y + mat2*y
    dy[end] = 2*(y[end-1] - y[end])

    dz[1] = 0
    dz[2:end-1] .= mat1*z + mat2*z
    dz[end] = 2*(z[end-1] - z[end])
end
u0 = zeros(N+2,2); u0[1,:] .= 1
tspan = (0.0,5.0)
prob = ODEProblem(f, u0, tspan)

# Accelerate it
sys = modelingtoolkitize(prob)
fastprob = ODEProblem(sys, u0, tspan, jac=true, sparse=true)

# Non-Stiff Choice
# 12.908 ms (1073 allocations: 9.58 MiB)
@btime solve(fastprob,Tsit5(),saveat=0.01,abstol=1e-10,reltol=1e-9)

# Unaccelerated Non-Stiff Choice
# 21.300 ms (9658 allocations: 36.36 MiB)
@btime solve(prob,Tsit5(),saveat=0.01,abstol=1e-10,reltol=1e-9)

# Stiff Choice
# 479.284 ms (21441 allocations: 619.75 MiB)
@btime solve(fastprob,TRBDF2(),saveat=0.01,abstol=1e-10,reltol=1e-9)

# Unaccelerated Stiff Choice
# 15.849 s (1381204 allocations: 385.17 MiB)
@btime solve(prob,TRBDF2(),saveat=0.01,abstol=1e-10,reltol=1e-9)

You can probably take that and modify it to your original problem.