How do MATLAB's ODE solvers use a functional form of the Jacobian?

How do the MATLAB ode*** solvers, in particular ode15s, use a functional form of the Jacobian? Does the functional form have to be an analytical form, meaning no real values are present in the Jacobian but are instead calculated inside the solver...OR...is the functional form of the Jacobian supposed to return specific values, i.e. a numerical solution? Furthermore, how often is the Jacobain calculated? Lastly, what are the expected input variables into the function supposed to be?

I have been looking for the past couple of weeks and the only examples I can find show the use of either a sparse matrix or a constant Jacobian. This includes the official page for odeset (<http://de.mathworks.com/help/matlab/ref/odeset.html)>.

I am simulating a system of 53 variables all of which have a dependence on multiples, but not all, of the variables. Since each of the variables have a time dependence, then each therefore each location of the Jacobian varies in time too.

I've created a function to calculate the full Jacobian at a given time step t. This function requires input of t and n(t), n being the solution variables; however, this slows down my simulations significantly compared to the Sparse Matrix.

Below are some small snip-its of my code. I don't intend to provide everything, because this is more a question of "How do I use...?" instead of "Why won't ... work?"

options1 = odeset(...,'Jacobian',@Discharge_Jacobian2,...);
...
...
sol = ode15s(@Discharge_ODE,[0 delTau],n_Init,options1);
...
...
function [J] = Discharge_Jacobian2(t,n)
...
J = zeros([NumSpec]);
for i = 1:NumSpec
  ...
  for j = 1:NumReac
    ...
    for k = 1:NumSpec
      if Stoich(j,k,1) ~= 0 || Stoich(j,k,2) ~= 0
        J(k,i) = J(k,i) + RR_i(j,i);
      end
    end
  end
end

Compared to, where S is attached.

options1 = odeset(...,'JPattern',S,...);
...
...
sol = ode15s(@Discharge_ODE,[0 delTau],n_Init,options1);

Answer 1

[Damjan Lasic]

You are correct - your analytical jacobian should be a function that will return the jacobian matrix after input of t and all solutions at that t.

I think in your case most likely your defined jacobian is incorrect. You can check this by pausing the script during ode15s calculation without your jacobian, and check the approximated jacobian in the solver workspace (it calculates it using odenumjac function). Then input the solution and time to your jacobian. If your jacobian function is correct, then you should get similar results of the jacobian as the ode15s computes it numerically. I had similar issues and always found out that the probem was in an incorrect provided jacobian.

Also, it's good to re-check the orientation of the resulting Jacobian matrix. Rows should be dF(1-n)/dx and columns should be dF/dx(1-n) if that makes sense. It's clearly defined in the odeset documentation. Another good way to double check if your jacobian is correct is to use it on a smaller system where you can write the correct jacobian by hand, and then compare the results of both jacobian function with some random input.

Another option is that your jacobian function evaluates the jacobian very slowly. In your case you are using multiple loops which may not be the fastest. Still, i think its most likely an incorrect jacobian and this should be the first thing to verify.

Kind regards, Damjan