Re: How to pass time-dependent vectors with lsode

[Juan Pablo Carbajal]

On Mon, Aug 8, 2016 at 8:02 PM, Ferdinando <address@hidden> wrote:
> Hi Everyone,
>
> I'm trying to solve this problem in Octave. I am working on a particle track
> data where I can extract the information of local intensity (I) and time.
> Now I want to pass these information to lsode where i calculate a first
> order concentration (C) decay like this: dCdt = -k*I*C.
> I have difficulties to pass the time-dependent vector I(t).
> Can anyone help me in this?
>
> Here an example:
>
> time=[0 1 2 3 4 5];
> I = [0 1 1 2 2 0];
> local_time = [0 1 1 1 1 1];
>
> C0 = 100;
>
> y=lsode(@myODE,C0,time,I);
>
>
> -------------------
> function dCdt=myODE(x,time,I)
> k = 0.01;
>
> C = x(1);
>
> dCdt = -k.*I.*C;
>
> endfunction
>
>
>
> --
> View this message in context: 
> http://octave.1599824.n4.nabble.com/How-to-pass-time-dependent-vectors-with-lsode-tp4679062.html
> Sent from the Octave - General mailing list archive at Nabble.com.
>
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://lists.gnu.org/mailman/listinfo/help-octave

Most good numerical integrators use internal time steps to give the
solution at the required timestamps. This means lsode needs values of
intensities at arbitrary time. One possibility is to pass a function
handle that interpolates the intensity.
(code is not checked)

function dCdt=myODE(x,time,If)
k = 0.01;
C = x(1);
dCdt = -k.*If(time).*C;
endfunction

time=[0 1 2 3 4 5];
I = [0 1 1 2 2 0];
pp = interp1 (time, I, "pp");
If =@(t) ppval (pp, t);

C0 = 100;
y=lsode(@myODE,C0,time,If)

Another option, is to use the analytical solution of your equation
when I is written as a piecewise constant function.

Cheers