How does MATLAB's ode15 handle integration of an equation containing a matrix of values?

Question

I am translating MATLAB code into Python, but before worrying about the translation I would like to understand how MATLAB and specifically its ODE15s solver are interpreting an equation.

I have a function script, which is called upon in the master script, and this function script contains the equation:

function testFun=testFunction(t,f,dmat,releasevec)

    testFun=(dmat*f)+(releasevec.');

Within testFunction, t refers to time, f to the value I am solving for, dmat to the matrix of constants I am curious about, and releasevec to a vector of additional constants.

The ODE15s solver in the master script works its magic with the following lines:

for i=1:1461

    [~,f]=ode15s(@(t, f) testFunction(t, f, ...
    [dAremoval(i), dFWtoA(i), dSWtoA(i), dStoA(i), dFSedtoA(i), dSSedtoA(i); ...
    dAtoFW(i), dFWremoval(i), dSWtoFW(i), dStoFW(i), dFSedtoFW(i), dSSedtoFW(i); ...
    dAtoSW(i), dFWtoSW(i), dSWremoval(i), dStoSW(i), dFSedtoSW(i), dSSedtoSW(i); ...
    dAtoS(i), dFWtoS(i), dSWtoS(i), dSremoval(i), dFSedtoS(i), dSSedtoS(i); ...
    dAtoFSed(i), dFWtoFSed(i), dSWtoFSed(i), dStoFSed(i), dFSedremoval(i), dSSedtoFSed(i); ...
    dAtoSSed(i), dFWtoSSed(i), dSWtoSSed(i), dStoSSed(i), dFSedtoSSed(i), dSSedremoval(i)], ...
    [Arelease(i), FWrelease(i), SWrelease(i), Srelease(i), FSedrelease(i), SSedrelease(i)]), [i, i+1], fresults(:, i),options);
    fresults(:, i + 1) = f(end, :).';

fresults is a table initially of zeros that houses the f results. The options call odeset to get 'nonnegative' values. The d values matrix above is a 6x6 matrix. I already have all of the d values and release value calculated. My question is: how is ode15s performing the integration with a 6x6 matrix given in the testfunction equation? I have tried to solve this by hand, but have not been successful. Any help would be much appreciated!!

#

def func(y, t, params):
    f = 5.75e-16
    f = y
    dmat, rvec = params
    derivs = [(dmat*f)+rvec]
    return derivs

# Parameters
dmat = np.array([[-1964977.10876756, 58831.976165, 39221.31744333,  1866923.81515922, 0.0, 0.0],
             [58831.976165, -1.89800738e+09, 0.0, 1234.12447489, 21088.06180415, 14058.70786944],
             [39221.31744333, 0.84352331, -7.59182852e+09, 0.0, 0.0, 0.0],
             [1866923.81515922, 0.0, 0.0, -9.30598884e+08, 0.0, 0.0],
             [0.0, 21088.10183616, 0.0, 0.0, -1.15427010e+09, 0.0],
             [0.0, 0.0, 14058.73455744, 0.0, 0.0, -5.98519566e+09]], np.float)
new_d = np.ndarray.flatten(dmat)

rvec = np.array([[0.0], [0.0], [0.0], [0.0], [0.0], [0.0]])

f = 5.75e-16

# Initial conditions for ODE
y0 = f

# Parameters for ODE 
params = [dmat, rvec]

# Times
tStop = 2.0
tStart = 0.0
tStep = 1.0

t = np.arange(tStart, tStop, tStep)

# Call the ODE Solver
soln = odeint(func, y0, t, args=(params,))
#y = odeint(lambda y, t: func(y,t,params), y0, t)

In a system with a 6x6 matrix and a 6x1 vector, the initial state y0 also needs to be a 6x1 vector. You are providing a scalar instead. — Lutz Lehmann

crbah crbah · Accepted Answer · 2015-07-10T06:31:04

It says here that ode15s uses backward difference formula for differentiation. Your differential equation is (as far as I understand) f' = testFunc(t,f) and it has some vector matrix calculations inside the function. Then you can replace the differentiation by a backward difference formula that is:

f_next = f_prev + h*testFunc(t,f_next);

where f_prev is the initial values of the vector. Here there is no important difference in calculations just because testFunc(t,f) function includes a 6x6 matrix. Each time it solves an inverse problem to find f_next by creating Jacobian matrices numerically.

However, trying to code algorithms as matlab does may be harder than we think since matlab has some (optimization related or not) special treatments to the problems. You should be careful on each value you get.

How does MATLAB's ode15 handle integration of an equation containing a matrix of values?

2 Answers