Update variables within an anonymous function? (Using ode45 on a discrete set of DE's)

Hello,

I am trying to use ode45 to implement 4th order Runge-Kutta on a system of first order differential equations. Within the body of my code, I call ode45 as follows:

    [times, T] = ode45(@(time, T) interpolate(time, A, ...
            int_data, tmax, temp_dist, a, b), ...
            subt, temp_dist);

where @(time, T) interpolate(...) is my function handle, subt is my interval of integration, and temp_dist is my initial condition. interpolate.m is defined as follows:

    function dTdt = interpolate(t, A, int_data, tmax, temp_dist, a, b)

    m = interp1([1:tmax], int_data, t, 'linear', 'extrap'); 
    % interpolate int_data at time t 

    m(m < 0) = 0;
    % none of the values in m can be negative

    % define matrix A at interpolated times (the first row of A depends on
    % t)
    A(1,1) = -1*b-a*m(3)-a*m(4);
    A(1,2) = b;
    A(1,end) = a*(m(1) + m(2) + m(3)*m(5) + ...
                    m(4)*m(5));

    % define output (evaluate system of ODE's at time t)
    dTdt = A*temp_dist;

where t is the time input by ode45, A is a matrix of coefficients whose first row changes with t, int_data is a matrix of information used to determine the first row of A, temp_dist is a vector of initial conditions, and tmax, a and b are used to determine the first row of A.

As you probably know, ode45 calls interpolate.m each time it needs a new value of dTdt for a specific time t. My problem is that temp_dist needs to be updated in between calls. After dTdt is calculated, temp_dist should be updated to:

    temp_dist = temp_dist + dTdt;

This new value of temp_dist should be the one used by ode45 the next time interpolate.m is called.

I'm stuck trying to implement this change. The obvious thing might be to redefine temp_dist within ode45.m, but the source code for this function is pretty incomprehensible to me. I've tried redefining temp_dist as shown above within interpolate.m, but when ode45 calls interpolate.m the next time, it's back to its original value.

Does anyone know how to solve this problem? My code is written in R2014a.

Thanks!!

[John D'Errico]

First of all, do NOT modify ODE45. Too much of a risk in that for you.

Next, it looks as if you are trying to do an update

temp_dist = temp_dist + dTdt;

that is essentially a simple linear step, as if you are trying to use ODE45, but making an Euler correction in the middle. This will simply degrade the performance of any 4th order method you employ, so a bad thing, and a waste of time to try to use ODE45.

Answer 1

[Torsten]

I guess that temp_dist are your solution variables. Thus your call to ODE45 has to be

[times, T] = ode45(@(time, temp_dist) interpolate(time, A, ...
            int_data, tmax, temp_dist, a, b), ...
            subt, temp_dist0);

where temp_dist0 are your initial conditions for temp_dist.

In this way, ODE45 will automatically update dempdist according to

temp_dist = temp_dist + dt*dTdt

But to be honest:

From what you did in your function "interpolate", I can neither identify what your solution variables for ODE45 are nor how your ODEs you are trying to solve look like.

Best wishes

Torsten.