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


i translating matlab code python, before worrying translation understand how matlab , ode15s solver interpreting equation.

i have function script, called upon in master script, , function script contains equation:

function testfun=testfunction(t,f,dmat,releasevec)      testfun=(dmat*f)+(releasevec.'); 

within testfunction, t refers time, f value solving for, dmat matrix of constants curious about, , releasevec vector of additional constants.

the ode15s solver in master script works magic 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(:, + 1) = f(end, :).'; 

fresults table of zeros houses f results. options call odeset 'nonnegative' values. d values matrix above 6x6 matrix. have of d values , release value calculated. question is: how ode15s performing integration 6x6 matrix given in testfunction equation? have tried solve hand, have not been successful. 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 ode y0 = f  # parameters ode  params = [dmat, rvec]  # times tstop = 2.0 tstart = 0.0 tstep = 1.0  t = np.arange(tstart, tstop, tstep)  # call ode solver soln = odeint(func, y0, t, args=(params,)) #y = odeint(lambda y, t: func(y,t,params), y0, t) 

it says here ode15s uses backward difference formula differentiation. differential equation (as far understand) f' = testfunc(t,f) , has vector matrix calculations inside function. can replace differentiation backward difference formula is:

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

where f_prev initial values of vector. here there no important difference in calculations because testfunc(t,f) function includes 6x6 matrix. each time solves inverse problem find f_next creating jacobian matrices numerically.

however, trying code algorithms matlab may harder think since matlab has (optimization related or not) special treatments problems. should careful on each value get.


Comments