Provide Jacobian Matrix for ODE15s solver

I'm using the ODE15s solver and want to ad the input for a jacobian matrix. The following file calculates the jacobi matrix J:

function [J] = Jacobi_Matrix(Fss,y_sym)
%JACOBI_MATRIX calculates the Jacobi Matrix for ODE15s Solver

%Create jacobi matrix

ny = length(y_sym);

y_sub = transpose(sym('y',[1 ny])); %substitute variables

Fss = subs(Fss,y_sym,y_sub);

J = jacobian(Fss,y_sub);
J = subs(J,{sym('u'), sym('v')}, {u, v});
J = subs(J, y_sub,y);

J = double(J)

end

Where y(1)...y(16) are the the degrees of freedom in state space representation and u,v are input data from earthquake excitation which are loaded in to Fss which is the right side of the Equation. The problem now is, that

   J = subs(J,{sym('u'), sym('v')}, {u, v});
   J = subs(J, y_sub,y);

is not working, it always says 'undefined variables'. u and v are loaded in my ODE file:

function dy = rocking_DGL(t,y,time_r,dt,ug,vg,Fss,Sys,y_sym)

r = Sys.r;
roh = Sys.roh;
A = Sys.A;
I = Sys.I;
ne = Sys.ne;
nnode = Sys.nnode;
ndof = nnode*3;
Le = Sys.Le;
E = Sys.E;
g = Sys.g;
EI = E*I;
EA = E*A;

%Interpolate Earthquake Excitation

i = floor(t/dt) + 1;
if i < length(time_r)
     u = ug(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(ug(i+1)-ug(i));
     v = vg(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(vg(i+1)-vg(i));

else
      u = 0;
      v = 0;

end

%Substitute Parameters

Fss = subs(Fss,y_sym,y);
Fss = subs(Fss,{sym('u'), sym('v')}, {u, v});
%Load State Space Representation of Problem

dy = double(Fss);

end

The substitution in the ODE function file for u and v works because they are defined there. What do I need to change to load the Jacobian into my odesolver options? What is the best command/way to tell the solver to use the jacobian?

Thanks for your answers!

[Marius]

I changed the code as followed:

refine = 1; 
    options = odeset('Refine', refine,...
        'RelTol', 1e-5,'absTol',1e-5*ones(1,length(y_sym)),'Jacobian',@JacobiFun,'JConstant','no','Mass',MassFun,'MassSingular','no','MstateDependence','strong','OutputFcn',@odetpbar);

      [t_up,y_up,te,ye,ie] = ode15s(@(t,y) ...
          rocking_DGL(t,y,time_r,dt,ug,vg,Fss,Sys,y_sym),time_r,IC_rocking,options);

Where the function JacobiFun is defined as:

function dfdy = JacobiFun(t,y,time_r,dt,ug,vg,J,y_sym,Sys)
%This Function evaluates the Jacobi Matrix

r = Sys.r;
roh = Sys.roh;
A = Sys.A;
I = Sys.I;
ne = Sys.ne;
nnode = Sys.nnode;
ndof = nnode*3;
Le = Sys.Le;
E = Sys.E;
g = Sys.g;
EI = E*I;
EA = E*A;

%Interpolate Earthquake Excitation

i = floor(t/dt) + 1;
if i < length(time_r)
     u = ug(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(ug(i+1)-ug(i));
     v = vg(i) + (t-time_r(i))/(time_r(i+1)-time_r(i))*(vg(i+1)-vg(i));

else
      u = 0;
      v = 0;

end

%Substitute Parameters

J = subs(J,y_sym,y);
J = subs(J,{sym('u'), sym('v')}, {u, v});

%Load State Space Representation of Problem

dfdy = double(J);

end

Where J is the Jacobi I already calculated before which is symbolic. When I try to start the solver it gives an Error:

Not enough input arguments.

Error in JacobiFun (line 4) r = Sys.r;

Error in ode15s (line 350) dfdy = feval(Jac,t,y,Jargs{:});

Why are there not enough arguments? I pass them with Sys!

Answer 1

[Torsten]

options=...,'Jacobian',@(t,y)JacobiFun(t,y,time_r,dt,ug,vg,J,y_sym,Sys),...

Best wishes

Torsten.

[Marius]

Thanks! It worked.