Setting up ode solver options to speed up compute time

Hi All,

I'm specifying the `'JPattern', sparsity_pattern` in the ode options to speed up the compute time of my actual system. I am sharing a

sample code below to show how I set up the system using a toy example. Specifying the `JPattern` helped me in reducing the compute time from 2 hours to 7 min for my real system. I'd like to know if there are options (in addition to `JPatthen`) that I can specify to further decrease the compute time . I found the `Jacobian` option but I am not sure how to compute the Jacobian easily for my real system.

global mat1 mat2
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;
        ];



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

f0 = fun(0, x0);
joptions = struct('diffvar', 2, 'vectvars', [], 'thresh', 1e-8, 'fac', []);
J = odenumjac(@fun,{0 x0}, f0, joptions); 
sparsity_pattern = sparse(J~=0.);
options = odeset('Stats', 'on', 'Vectorized', 'on', 'JPattern', sparsity_pattern);

ttic = tic();
[t, sol]  =  ode15s(@(t,x) fun(t,x), tspan , x0, options);
ttoc = toc(ttic)
fprintf('runtime %f seconds ...\n', ttoc)
plot(t, sol)

function f = fun(t,x)
global mat1 mat2
% f = zeros('like', x)
% size(f)
f = zeros(size(x), 'like', x);
size(f);
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

Are there inbuilt options available for computing the Jacobian?

I tried something like the below

x = sym('x', [5 1]);
s = mat1*x + mat2*x;
J1 = jacobian(s, x)

But this takes huge time for large system.

Suggestions will be really appreciated.

Side note:

I would also like to know if there is someone on the forum to whom I can demonstrate my code and seek help to resolve the issue mentioned above.

Unfortunately, I cannot post my actual system here .

[Deepa Maheshvare]

Any suggestions?

[Bjorn Gustavsson]

When you write "Unfortunately, I cannot post my actual system here" you shout out most interested helpers. That phrase means that one would try to problem-shoot irrelevant toy-model-problems that will not help solving your real problem. That entire process will take much much more time to solve than necessary - this is frustrating for the kind souls that volunteer. You will now have to show enough interest in getting help by constructing a similar enough problem from a numerical standpoint that you share (same level of stiffness same everything else too except your "secret" parts.).

[Torsten]

In your toy example, the Jacobian is constant over time. Is this also true for your real system ?

[Deepa Maheshvare]

@Torsten Yes, it is true for my real system as well.

[Torsten]

Then calculate the Jacobian J once before calling ODE15S and pass it to the solver by

odeset('Jacobian',J)

Calculation of the Jacobian can be done using jacobianest from the file exchange (or using odenumjac as you did).

[Deepa Maheshvare]

Hi @Torsten, thank you. Could you please let me know if it is recommended to specify both `JPattern` and `Jacobian` in odeset?

[Torsten]

My guess is that specifying the pattern additionally to specifying the Jacobian itself might even slow down the solver because the nonzero positions have to be addressed in a complicated way. At least the gain in performance will be low in my opinion.

But what does your test example indicate ? I'm curious ...

[Deepa Maheshvare]

@Torsten, For my real system,

jpattern: 230 s

jpattern and jacobian: > 900 s

jacobian: > 600 s (still running)

[Torsten]

This can't be true in my opinion. jpattern + jacobian must be much more efficient than jpattern alone.

Or it was not true that the jacobian remains constant over time, as you claimed.

What about the results for the toy example ?

[Deepa Maheshvare]

@Torsten, please find below the runtimes for the toy exmaple. This works as expected i.e. specifying the jacobian gives better performance over `jpattern` or `jacobian+jpattern`. I'll try figuring out what's going on with my actual problem.

    global advMat diffMat
    
    gridSize = 10000;
    x0 = [1; zeros(gridSize - 1, 1)];
    tspan = 0:0.01:5;
    
    % Specify the finite difference matrices including boundary conditions
    % You'll probably need to add some 1/dx or 1/dx^2 scaling to the matrices
    e = ones(gridSize, 1);
    advMat = spdiags([e, -e], 0:1, gridSize, gridSize);
    diffMat = spdiags([e, -2*e, e], -1:1, gridSize, gridSize);
    
    f0 = f(0, x0);
    joptions = struct('diffvar', 2, 'vectvars', 2, 'thresh', 1e-8, 'fac', []);
    J = odenumjac(@f,{0 x0}, f0, joptions); 
    jpattern = sparse(J~=0.);
    
    % The problem is linear so the Jacobian is the sum of linear operators in RHS
    jacobian = advMat + diffMat;
    
    
    ttic = tic();
    % [t, sol] = ode15s(@(t, y) f(t, y), tspan, x0, odeset('Jacobian', jacobian));
    % [t, sol] = ode15s(@(t, y) f(t, y), tspan, x0, odeset('JPattern', jpattern));
    [t, sol] = ode15s(@(t, y) f(t, y), tspan, x0, odeset('Jacobian', jacobian, 'JPattern', jpattern));
    
    ttoc = toc(ttic);
    fprintf('runtime %f seconds ...\n', ttoc);
    plot(sol);
    
    function dy = f(~, y)
    global advMat diffMat
    dy = advMat * y + diffMat * y;
    end

Result:

For a grid size of 10000

    Jacobian: runtime 15.508172 seconds ...
    Jpattern: runtime 57.470325 seconds ...
    Jacobian + Jpattern: runtime 30.028399 seconds ...

For a grid size of 5000

    Jacobian: runtime 0.203265 seconds ...
    Jpattern: runtime 0.650601 seconds ...
    Jacobian + Jpattern: runtime 0.256899 seconds ...

For a grid size of 50000

    Out of memory. Type "help memory" for your options.

:/

[Torsten]

I'm quite certain that the Jacobian changes in time.

If this is the case, the Jacobian supplied becomes more and more inexact amd increasingly smaller time steps are required to compensate the error in the Jacobian.

[Deepa Maheshvare]

@Torsten

In my real system I am solving the advection diffusion dynamics. The matrices `mat1` and `mat2` presented in the original post are the advection and diffusion operators. Since these are constant matrices, the jacobian is a constant matrix in my real problem too. Please let me know if you sense this is wrong

.

[Torsten]

I don't know about your problem and implementation. You could test by making random vectors for your solution variables and make Matlab generate the corresponding Jacobians.

[Deepa Maheshvare]

Thank you

[Deepa Maheshvare]

Hi @Torsten,

"You could test by making random vectors for your solution variables and make Matlab generate the corresponding Jacobians."

I tried this for my toy example. Could you please have a look at the code posted below?

isequal(J, J1)

returns 0. I am not sure why J is not equal to J1.

global advMat diffMat

gridSize = 10;
x0 = [1; zeros(gridSize - 1, 1)];
tspan = 0:0.01:5;

e = ones(gridSize, 1);
advMat = spdiags([e, -e], 0:1, gridSize, gridSize);
diffMat = spdiags([e, -2*e, e], -1:1, gridSize, gridSize);

f0 = f(0, x0);
joptions = struct('diffvar', 2, 'vectvars', 2, 'thresh', 1e-8, 'fac', []);
J = odenumjac(@f,{0 x0}, f0, joptions)
jpattern = sparse(J~=0.);

x1 = x0+ rand(length(x0),1);
f1 = f(0, x1);
J1 = odenumjac(@f,{0 x1}, f1, joptions)
isequal(J, J1)

ttic = tic();
[t, sol] = ode15s(@(t, y) f(t, y), tspan, x0, odeset('Jacobian', J));
% [t, sol] = ode15s(@(t, y) f(t, y), tspan, x0, odeset('JPattern', jpattern));
% [t, sol] = ode15s(@(t, y) f(t, y), tspan, x0, odeset('Jacobian', jacobian, 'JPattern', jpattern));

ttoc = toc(ttic);
fprintf('runtime %f seconds ...\n', ttoc);
plot(sol);

function dy = f(~, y)
global advMat diffMat
dy = advMat * y + diffMat * y;
end

[Torsten]

de.mathworks.com/matlabcentral/fileexchange/15816-isalmost

[Deepa Maheshvare]

@Torsten, Thank you, this function is a really useful one. I checked using

r = isalmost  (J, J1, 1e-3)

Still all entries of `r` isn't 1. The entries (1,2) and (2,2) are different. Unfortunately, I couldn't understand why this difference occurs for a constant jacobian matrix.

[Torsten]

Analytical Jacobian should be Jac_ana = advMat + diffMat.

Maybe you can just output J, J1 and Jac_ana and compare them directly.

[Bjorn Gustavsson]

@Deepa Maheshvare - if you're solving a diffusion-advection problem then maybe it is worthwhile to look at the PDE-solvers, if you have access to the pde-toolbox.