help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: How to pass time-dependent vectors with lsode


From: Juan Pablo Carbajal
Subject: Re: How to pass time-dependent vectors with lsode
Date: Mon, 8 Aug 2016 23:25:24 +0200

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



reply via email to

[Prev in Thread] Current Thread [Next in Thread]