Setting up ode solver options to speed up compute time

8 次查看（过去 30 天）

Deepa Maheshvare 2021-5-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/833253-setting-up-ode-solver-options-to-speed-up-compute-time

评论： Bjorn Gustavsson 2021-5-25

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 .

19 个评论
显示 17更早的评论隐藏 17更早的评论

Deepa Maheshvare 2021-5-22

在 MATLAB Online 中打开

@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.

Deepa Maheshvare 2021-5-23

在 MATLAB Online 中打开

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 2021-5-23

Analytical Jacobian should be Jac_ana = advMat + diffMat.

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

Bjorn Gustavsson 2021-5-25

@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.

请先登录，再进行评论。

请先登录，再回答此问题。