TransWikia.com

Solving differential equation by specifying jacobian pattern

Computational Science Asked on February 28, 2021

This is a follow up to my previous question posted here

I’m trying to construct the sparsity pattern of the jacobian matrix to speed up the computation of a large system of odes. The following is the code in which I am trying to set up the jpattern in
odeset for a toy model in MATLAB.

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', 1, 'thresh', 1e-8, 'fac', []);
J = odenumjac(@fun,{0 x0}, f0, joptions); 
sparsity_pattern = sparse(J~=0.);
options = odeset('Stats', '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

fprintf('size of x: %d %dn', size(x))

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

When I run the code, I notice that the size of the vector x in fun changes for the second iteration
because of which an error occurs.

size of x: 10 1
size of x: 10 10

Error:

Unable to perform assignment because the size of the left side is 8-by-1 and the size of the right side is 8-by-10.

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

Error in odenumjac (line 143)
    Fdel = feval(F,Fargs_expanded{:});

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

Suggestions on how to fix this error will be really helpful.

EDIT:

I’ve tried to vectorize f in fun and also set vectvars=2 to vectorize the jacobian calculation.

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', 2, 'thresh', 1e-8, 'fac', []);
J = odenumjac(@fun,{0 x0}, f0, joptions); 
sparsity_pattern = sparse(J~=0.);
options = odeset('Stats', 'on', 'JPattern', sparsity_pattern, 'Vectorized', 'on');

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(1,:) = 0;
f(2:9,:) = mat1*x + mat2*x;
f(10,:) = 2*(x(end-1) - x(end));
end

However, there is a problem again when vectvars = 2 in joptions and/or ‘Vectorized’, ‘on’ in options defined for ode15s.

Unable to perform assignment because the size of the left side is 8-by-1 and the size of the right side is 8-by-10.

Error in cse_11_5_20>fun (line 44)
f(2:9,:) = mat1*x + mat2*x;

Error in odenumjac (line 143)
    Fdel = feval(F,Fargs_expanded{:});

Error in cse_11_5_20 (line 31)
J = odenumjac(@fun,{0 x0}, f0, joptions);

One Answer

odenumjac calls your function in a vectorized manner it seems, and your function is not vectorized. You can easily change that by changing the second index of f in your function to : instead of 1, for instance: f(10,:) = 2*(x(end-1,:) - x(end,:));

I thought the setting joptions.vectvars=1 would not allow the vectorised call (see one of your other questions). I realize that this is actually not the case, you should instead set it to [] (empty). If you want to vectorize the Jacobian calculation, set it to 2. You can type open odenumjacin the Matlab console to access the function file and read its documentation.

More info on vectorization can be found here: https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html

Correct answer by Laurent90 on February 28, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP