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

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.

• Provide the analytical Jacobian (or at least a sparsity pattern). Is your ODE stiff? If not, try an explicit method (e.g., ode45). How do Matrix1 and Matrix2 look like? Do Y and Z interact (if not, solve two smaller ODEs). – cos_theta Nov 15 '20 at 23:38
• Can you share the MATLAB code? I'd like to demonstrate the speed in Julia, which is generally between 20x-100x on this kind of thing. Of course, it would be hard to show in your case without seeing the code. – Chris Rackauckas Nov 28 '20 at 8:22
• @BillGreene Could you please check my edits? I couldn't share my actual code here (which has convection + diffusion + reaction). I tried setting up a toy model and the code snippet is shared in the original post – Natasha Feb 16 at 17:39
• @ChrisRackauckas Could you please check my edit? – Natasha Feb 16 at 17:40

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

• Regarding vectorisation, I read in the only answer to this question in the MATLAB forum, that "it is a bad idea from a numerical point of view [...]. When you run the integration over a set of different initial conditions, the step size is limited by the most sensitive coordinate and therefore the accumulated discretization error is much higher than possible" and he proceeds to further justification. Do you think for that reason we should aim to avoid it? – kostas1335 Nov 26 '20 at 19:35
• For me it is not clear from the original post that the OP wants to run multiple simulations in parallel. In that case, yes, the time step will be limited by the "hardest" solution... An alternative is to only run one simulation at a time, but using different threads to compute multiple solutions in parallel. I'll edit my answer with that. – Laurent90 Nov 26 '20 at 19:48
• @Laurent90 sorry for the really late reply. Thanks a lot for the detailed explanation. Yes, my problem is stiff, variables Y and Z are coupled and the total number of variables is around 2400. I could change Matrix1 and Matrix2 as sparse matrices. Unfortunately, it's not clear to me how the Jpattern has to be specified. Could you please direct me to some sample problems from where I can learn how to set this up? – Natasha Feb 11 at 7:08
• @Natasha In your case, you can easily find the form of the Jacobian of your system (shoudl be something like a block-diagonal matrix, with two equal blocks, which are diag(vector_1)*(matrix1+matrix2). In a more general case, you can compute your Jacobian via finite-differences (see the function odenumjac in Matlab for instance). I sometimes do this once to compute the sparsity pattern of my Jacobian, and then give that Jacobian pattern to the Matlab ODE solver, so that it is able to more efficiently compute the Jacobian when it needs to. – Laurent90 Feb 12 at 17:06
• @Natasha I have updated my answer with more information regarding what you ask. I think you should look up more in-depth explanations if you want to have more precise knowledge ;) – Laurent90 Feb 15 at 18:07

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.

• Thanks a lot! I'd love to know if it is possible to call Julia from MATLAB. For instance, I'm interested to know if the MATLAB function fun in my code can be passed as an argument here sys = modelingtoolkitize(prob)` . – Natasha Feb 17 at 7:58
• There are ways to do it but they aren't well-maintained. I wouldn't recommend trying it except for fun. Also, it's going to hurt performance if you do it like that. – Chris Rackauckas Feb 17 at 15:30