Why am I not able to run this code?

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

I've defined the `options` for `odenumjac` and it still throws the error

Not enough input arguments.

Am I missing something?

Answer 1

[KALYAN ACHARJYA]

Here

J = odenumjac(fun,{0 x0}, f0, joptions); 
^..............^ function without input arguments 

[Deepa Maheshvare]

Could you please explain your answer a bit? Unfortunately, it is not clear to me what's wrong in my code (input arguments to `fun` is missing?) and how the function `F` has to be passed to `odenumjac(F,Fargs,Fvalue,options)`.

In the description given for `odenumjac` it is specified that

The arguments of F are specified in a cell array FARGS

So I've tried

fun,{0 x0}

[Walter Roberson]

J = odenumjac(@fun,{0 x0}, f0, joptions);

[Deepa Maheshvare]

Hi @Walter Roberson,

I gave this `J = odenumjac(@fun,{0 x0}, f0, joptions);` a try .

Unfortunately, the function handle didn't help

Undefined function 'odenumjac' for input arguments of type 'function_handle'.

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

[Walter Roberson]

toolbox/matlab/funfun/private/odenumjac.m  % Private to funfun

public use of odenumjac seem to disappear after 2008. The function still exists but it is private now.

[KALYAN ACHARJYA]

Is that function "odenumjac" MATLAB inbuilt ?

[Deepa Maheshvare]

@KALYAN ACHARJYA

I did

>> exist odenumjac

ans =

     0

I am not sure if this means that the function `odenumjac` is not in-built

[Deepa Maheshvare]

@Walter Roberson

I'll happy to know if there is a workaround for this. I'm trying to use odenumjac to find the sparsity pattern of jacobian. My system is pretty large with ~9000 odes modelling convection diffusion reaction . Without the jacobian , it takes ~2 hours to simulate my model.

[Walter Roberson]

which -all odenumjac

You will see it in a /private directory. Copy it to your working directory.

You might need other functions from the same directory.

[Deepa Maheshvare]

@Walter Roberson

Thanks so much. That helps!

Unfortunately, I still couldn't compute the sparsity pattern successfullly.

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 %d\n', size(x))

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

The above is the updated code and following issue occurs

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);
 

I see that the size of `x` changes in the second iteration.

>> test
size of x: 10 1
size of x: 10 10

I'm not sure why this happens. Suggestions on how this issue can be fixed will be of great help.

[KALYAN ACHARJYA]

For this issue

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

Once I tried to reproduce the error with a little modification to the function file, I could not find any such issue. Note that here I only partially executed the main code upto the following line

f0 = fun(0, x0);

Function file: Input t is not used within the function

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

[Walter Roberson]

joptions = struct('diffvar', 2, 'vectvars', 1, 'thresh', 1e-8, 'fac', []);

That 'vectvars', 1 option says that your function is vectorized over the first variable, t,

%   ODENUMJAC takes advantage of vectorization, i.e., when several values F 
%   can be obtained with one function evaluation. Set OPTIONS.VECTVAR 
%   to the indices of vectorized arguments: VECTVAR = [2] indicates that 
%   F(t,[x1 y2 ...]) returns [F(t,x1) F(t,x2) ...], while VECTVAR = [1,2] 
%   indicates that F([t1 t2 ...],[x1 x2 ...]) returns [F(t1,x1) F(t2,x2) ...]. 

So you have declared that it an pass in multiple t values and the function will return a result for each of the t values -- but your code does not in fact handle that situation.