octave and delay ode's

[franco basaglia]

Some parts of this message have been removed. Learn more about Nabble's security policy.

Hi all,

I'm trying to solve a system of Seven ODE's (from Stella) by Octave using lsode.
Two of this are delay differential equations. For example on state variable S(t):

%S(t)
% S(t)=S(t−dt)+(dSi−dSo)×dt

% dSi = 0.2×DELAY((F+0.5×U)×(1−GP),5)×(1−S)×(1−T)×AE %this is Stella definition
% DELAY on S(t)
%%% F=x(2); U=x(4); GP=1-exp(d*x(5)); S=x(6); T=x(7)

if (t < 5*delta_t)
t_back = t;
t = 0;
delay_si = ((x(2)+0.5*x(4))*(exp(d*x(5))));
t = t_back;
t;
else
t_back = t;
t = t-5*delta_t;
delay_si = ((x(2)+0.5*x(4))*(exp(d*x(5))));
t = t_back;
endif

dSi = 0.2*delay_si*(1-x(6))*(1-x(7))*AE;

dSo = (0.1 + 1 - exp(d*x(5))) * (0.1 + x(7)) * x(6) * 0.2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xdot(6) = dSi - dSo ;


Using lsode(adaptive step size method), for Cycle doesn't work because there is an “incompatibility” on delta_t.
Maybe if I print a matrix with x(t) value for every dt I can run for Cycle on these values.
But I don't know how to print this matrix with lsode.

Any idea ?
Alternative way to solve ode's delay?

Thanks

best regards
f.b.

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

On 28 Oct 2009, at 11:07, franco basaglia wrote:

> Hi all,
>
> I'm trying to solve a system of Seven ODE's (from Stella) by Octave  
> using lsode.
> Two of this are delay differential equations. For example on state  
> variable S(t):
>
> %S(t)
> % S(t)=S(t−dt)+(dSi−dSo)×dt
>
> % dSi = 0.2×DELAY((F+0.5×U)×(1−GP),5)×(1−S)×(1−T)×AE  
> %this is Stella definition
> % DELAY on S(t)
> %%% F=x(2); U=x(4); GP=1-exp(d*x(5)); S=x(6); T=x(7)
>
> if (t < 5*delta_t)
> t_back = t;
> t = 0;
> delay_si = ((x(2)+0.5*x(4))*(exp(d*x(5))));
> t = t_back;
> t;
> else
> t_back = t;
> t = t-5*delta_t;
> delay_si = ((x(2)+0.5*x(4))*(exp(d*x(5))));
> t = t_back;
> endif
>
> dSi = 0.2*delay_si*(1-x(6))*(1-x(7))*AE;
>
> dSo = (0.1 + 1 - exp(d*x(5))) * (0.1 + x(7)) * x(6) * 0.2;
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> xdot(6) = dSi - dSo ;
>
>
> Using lsode(adaptive step size method), for Cycle doesn't work  
> because there is an “incompatibility” on delta_t.
> Maybe if I print a matrix with x(t) value for every dt I can run for  
> Cycle on these values.
> But I don't know how to print this matrix with lsode.
>
> Any idea ?
> Alternative way to solve ode's delay?
>
> Thanks
>
> best regards
> f.b.

I don't really understand your code but for solving DDEs you might  
want to take a look at the ode**d [1] in the octave-forge package  
"odepkg" [2].
c.

[1] http://octave.sourceforge.net/doc/f/ode23d.html
     http://octave.sourceforge.net/doc/f/ode45d.html
     http://octave.sourceforge.net/doc/f/ode54d.html
     http://octave.sourceforge.net/doc/f/ode78d.html

[2] http://octave.sourceforge.net/odepkg/index.html
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[franco basaglia]

""I don't really understand your code""

Pratically I have to convert to Octave this Stella equation:
 dSi = 0.2×DELAY((F+0.5×U)×(1−GP),5)×(1−S)×(1−T)×AE

where (F+0.5×U)×(1−GP) lag behind by 5 time units. For the first 5 time units of the simulation, the delay will return the initial values.

So I use for cycle.


""for solving DDEs you might want to take a look at the ode**d [1] in the octave-forge package "odepkg" [2]""

I know this package but I don't know how to apply to my system. I'm a newbie. Have you some examples?


thanks

best regards
f.b.

2009/10/28 Carlo de Falco <carlo.defalco@...>

On 28 Oct 2009, at 11:07, franco basaglia wrote:

Hi all,

I'm trying to solve a system of Seven ODE's (from Stella) by Octave using lsode.
Two of this are delay differential equations. For example on state variable S(t):

%S(t)
% S(t)=S(t−dt)+(dSi−dSo)×dt

% dSi = 0.2×DELAY((F+0.5×U)×(1−GP),5)×(1−S)×(1−T)×AE %this is Stella definition
% DELAY on S(t)
%%% F=x(2); U=x(4); GP=1-exp(d*x(5)); S=x(6); T=x(7)

if (t < 5*delta_t)
t_back = t;
t = 0;
delay_si = ((x(2)+0.5*x(4))*(exp(d*x(5))));
t = t_back;
t;
else
t_back = t;
t = t-5*delta_t;
delay_si = ((x(2)+0.5*x(4))*(exp(d*x(5))));
t = t_back;
endif

dSi = 0.2*delay_si*(1-x(6))*(1-x(7))*AE;

dSo = (0.1 + 1 - exp(d*x(5))) * (0.1 + x(7)) * x(6) * 0.2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xdot(6) = dSi - dSo ;


Using lsode(adaptive step size method), for Cycle doesn't work because there is an “incompatibility” on delta_t.
Maybe if I print a matrix with x(t) value for every dt I can run for Cycle on these values.
But I don't know how to print this matrix with lsode.

Any idea ?
Alternative way to solve ode's delay?

Thanks

best regards
f.b.

I don't really understand your code but for solving DDEs you might want to take a look at the ode**d [1] in the octave-forge package "odepkg" [2].
c.

[1] http://octave.sourceforge.net/doc/f/ode23d.html
   http://octave.sourceforge.net/doc/f/ode45d.html
   http://octave.sourceforge.net/doc/f/ode54d.html
   http://octave.sourceforge.net/doc/f/ode78d.html

[2] http://octave.sourceforge.net/odepkg/index.html

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

On 28 Oct 2009, at 12:02, franco basaglia wrote:

> ""I don't really understand your code""
>
> Pratically I have to convert to Octave this Stella equation:
>  dSi = 0.2×DELAY((F+0.5×U)×(1−GP),5)×(1−S)×(1−T)×AE
>
> where (F+0.5×U)×(1−GP) lag behind by 5 time units. For the first  
> 5 time units of the simulation, the delay will return the initial  
> values.
> So I use for cycle.

sorry, I don't know what "stella" is so I do not understand this syntax,
maybe you could describe your system in plain mathematical syntax?

>
> ""for solving DDEs you might want to take a look at the ode**d [1]  
> in the octave-forge package "odepkg" [2]""
>
> I know this package but I don't know how to apply to my system. I'm  
> a newbie. Have you some examples?

yes, odepkg comes with a very nice set of examples and a pdf manual,  
type

odepkg_examples_dde ()

to see some DDE examples.
If you have installed odepkg, the manual should be located in the folder

~/octave/odepkg-0.6.9/doc

> thanks
>
> best regards
> f.b.

HTH,
c.
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[franco basaglia]

sorry, I don't know what "stella" is so I do not understand this syntax,

Stella is a simulation development environment. Delay sintax is in [1]


maybe you could describe your system in plain mathematical syntax?

I attach a pdf with my ode's system. Only equations number 6 and number 7 have a delay as you can see.




yes, odepkg comes with a very nice set of examples and a pdf manual, type
odepkg_examples_dde ()
to see some DDE examples.


Thanks.I'm trying to use it. I dont' yet understand how to apply in a system with
DDE's and ODE's



HTH,
c.

thanks for your collaboration

f.t.


[1]http://www.iseesystems.com/community/support/documentation/builtins/discrete_functions.htm#delay

2009/10/28 Carlo de Falco <carlo.defalco@...>

On 28 Oct 2009, at 12:02, franco basaglia wrote:

""I don't really understand your code""

Pratically I have to convert to Octave this Stella equation:
 dSi = 0.2×DELAY((F+0.5×U)×(1−GP),5)×(1−S)×(1−T)×AE

where (F+0.5×U)×(1−GP) lag behind by 5 time units. For the first 5 time units of the simulation, the delay will return the initial values.
So I use for cycle.

sorry, I don't know what "stella" is so I do not understand this syntax,
maybe you could describe your system in plain mathematical syntax?

""for solving DDEs you might want to take a look at the ode**d [1] in the octave-forge package "odepkg" [2]""

I know this package but I don't know how to apply to my system. I'm a newbie. Have you some examples?

yes, odepkg comes with a very nice set of examples and a pdf manual, type

odepkg_examples_dde ()

to see some DDE examples.
If you have installed odepkg, the manual should be located in the folder

~/octave/odepkg-0.6.9/doc

thanks

best regards
f.b.

HTH,
c.

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[franco basaglia]

I must be missing something as in your file I see no time derivatives...

sorry for my oversight. I attach again pdf with right system.

function f = fun (t, y, yd)
f(1) =-y(1) + yd(2);
f(2) =-y(2) ;
endfunction

ode23d (@fun, [0 5], [1;1], .1, [1;1])

notice that in the example above only one equation includes a delay:
is this what you mean by system of ODEs and DDEs?

Yes,sure.
To be more precise I have a system like this:

d y1(t)/ dt = -y1(t)

d y2(t)/ dt = -y2(t) + y1(t-5)
d y3(t)/dt  = -y3(t) + y2(t-10)*y1(t-10)

that I try to solve in this way:

function f = fun (t, y, yd)
f(1) =-y(1)

f(2) =-y(2) + yd(1)
f(3) =-y(3) + yd(2)*yd(1)
endfunction
t = [0:.5:20]
res = ode45d (@fun, t, [1;1;1], [5;10], ones (3,10))

It's run. But I don't think this is the right way.I don't know if with this notation it applies a delay of 5 time units to f(2) and of 10 time units to f(3). Note that in f(3) yd(1) comes again but with a different delay.
Infact for my script, that I attach here, is the same. It's runs but with different results compare to Stella.

Any Idea?

thanks a lot

f.t.

2009/10/29 Carlo de Falco <carlo.defalco@...>

On 29 Oct 2009, at 16:38, franco basaglia wrote:

sorry, I don't know what "stella" is so I do not understand this syntax,

Stella is a simulation development environment. Delay sintax is in [1]


maybe you could describe your system in plain mathematical syntax?

I attach a pdf with my ode's system. Only equations number 6 and number 7 have a delay as you can see.

I must be missing something as in your file I see no time derivatives...

yes, odepkg comes with a very nice set of examples and a pdf manual, type
odepkg_examples_dde ()
to see some DDE examples.


Thanks.I'm trying to use it. I dont' yet understand how to apply in a system with
DDE's and ODE's

ode23d (fun, [t0, tend], y0, tau, h)
solves the DDE sysytem

dy/dt = fun (t, y(t), y(t-tau))
y(t0) = y0
y(t) = h, t<t0

where of course y,y0,h can be vectors.

for example to solve the following system in 0<t<=5:

d y1(t)/ dt = - y1 (t) + y2 (t-.1)
d y2(t)/ dt = - y2 (t)

y1(t) = y2(t) = 1, t<=0

you can do:

function f = fun (t, y, yd)
f(1) =-y(1) + yd(2);
f(2) =-y(2) ;
endfunction

ode23d (@fun, [0 5], [1;1], .1, [1;1])

notice that in the example above only one equation includes a delay:
is this what you mean by system of ODEs and DDEs?

HTH,
c.

[scriptode45.m]

---

%%%%valori iniziali delle coperture
%copertura iniziale di prato semplice(tra 0 e 1)
global init_l = 0.148;
%copertura iniziale di prato incolto(tra 0 e 1)
global init_f = 0.147;
%copertura iniziale di prato complesso(tra 0 e 1)
global init_m = 0.303;
%copertura iniziale di arbusti(tra 0 e 1)
global init_s = 0.123;
%copertura iniziale di alberi(tra 0 e 1)
global init_t = 0.402;
%copertura iniziale di sottobosco(tra 0 e 1)
global init_u = (1 - init_l - init_f - init_m);

%%%%variabili di controllo
%altitudine (m)
global A = 600;
%Carico Animale Globale (UBA giorni/ha/anno)
global  GSD =103;

%%%%parametri
%Valore Pastorale medio incolto
global PVF = 10;
%Valore Pastorale medio pascolo semplice
global PVL = 20;
%Valore Pastorale medio pascolo complesso
global PVM = 40;
%Valore Pastorale medio sottobosco
global PVU = 5;
%Valore Pastorale globale
global GPV = 21;

%%%%equazioni
%valore pastorale locale
global LPV = PVF*init_f+PVL*init_l+PVM*init_m+PVU*init_u;
%carico animale locale(UBA giorno/ha/anno)
global init_lsd = GSD*LPV/GPV ;
%parametri
global d = -0.013305;
global b = -0.005;
global c = 1081.093;
%effetti dell'altitudine(tra 0 e 1)
global AE = exp(b*(A-c))/(exp(b*(A-c))+1);
%capacità portante globale(UBA giorno/ha/anno)
global GCC = 115000*GPV/18/A;
%tasso di utilizzazione globale(tra 0 e 1)
global GU =GSD/GCC;


%%%%%%%% parametri temporali di simulazione
%t inizio simulazione
global t_initial = 0;
%t fine simulazione
global t_final = 200;
       
global monitor_points = 400;
global delta_t = (t_final - t_initial) / monitor_points;


function f = fun (t,x,xd)

global d
global GSD
global PVF
global PVL
global PVM
global PVU
global GPV
global b
global A
global c
global AE
global SCr
global SCs
global SCp
global TCr
global TCs
global TCp
global GCC
global GU
global init_t
global t_initial
global t_final
global monitor_points
global delta_t

%x1
global init_l  
%x2
global init_f
%x3
global init_m
%x4
global init_u
%x5
global LPV
%x6
global init_lsd
%x7
global init_s

%%%%%% L(t) = L(t−dt)+(−LtoM−LtoF)×dt Prato Semplice
%%%%%% L(t)= -((0.06*L*GP)-(0,05*M*(1-GP))) - ((0,5*L*(1-GP)*(O.3+S)*AE)-(0.1*F*GP))
f(1) = -((0.06*x(1)*(1-exp(d*x(5))))-(0.05*x(3)*(exp(d*x(5)))))-( (0.5*x(1)*(exp(d*x(5)))*(0.3+x(6)))*AE-(0.1*x(2)*(1-exp(d*x(5)))));

%%%%%% F(t) = F(t−dt)+(LtoF+MtoF−FtoU)×dt  Incolto
%%%%%% F(t)= ((0,5*L*(1-GP)*(O.3+S)*AE)-(0.1*F*GP)) + (O.8*M*(1-GP)*(O.3+S)-(0.4*F*(GP)^2) - (F*(T-U))
f(2) =(((0.5*x(1)*(exp(d*x(5)))*(0.3+x(6)))*AE-(0.1*x(2)*(1-exp(d*x(5))))) + (0.8*x(3)*(exp(d*x(5)))*(0.3+x(6))-(0.4*x(2)*((1-exp(d*x(5))^2)))) - (x(2)*(x(6)-x(4))));

%%%%%% M(t) = M(t−dt)+(LtoM−MtoF)×dt Prato complesso
%%%%%% M(t)=((0.06*L*GP)-(0,05*M*(1-GP))) - (O.8*M*(1-GP)*(O.3+S)-(0.4*F*(GP)^2)
f(3) = ((0.06*x(1)*(1-exp(d*x(5))))-(0.05*x(3)*(+exp(d*x(5))))) - (0.8*x(3)*(exp(d*x(5)))*(0.3+x(6))-(0.4*x(2)*((1-exp(d*x(5))^2))));

%%%%%% U(t) = U(t−dt)+(FtoU)×dt
%%%%%% U(t) = (F*(T-U))
f(4) = (x(2)*(x(6)-x(4)));

%%%%% if per calcolare LSD
LCC =LPV*115000/18/A;
LU = x(5)/LCC;
if (LU < 1)
  dLSD = x(5) * (GU-LU);
else
  dLSD = x(5) *(1-LU);
endif

%%%%%% LSD(t) =LSD(t−dt)+(dLSD)×dt capacità di carico locale
f(5) = dLSD ;



%%%%% S(t)=S(t−dt)+(dSi−dSo)×dt        Arbusti

dSo = (0.1 + 1 - exp(d*x(5))) * (0.1 + x(7)) * x(6) * 0.2;
%
%%%%% funzione DELAY sulla variabile S(t)
%%%%% dSi = 0.2×DELAY((F+0.5×U)×(1−GP),5)×(1−S)×(1−T)×AE -----definizione di Stella
%
dSi = 0.2*((xd(2)+0.5*xd(4))*(exp(d*xd(5))))*(1-x(6))*(1-x(7))*AE;


f(6) = dSi - dSo;



%%%%% T(t) = T(t−dt)+(dTi−dTo)×dt   Alberi

dTo = x(7)*0.005;
%
%%%%% funzione DELAY sulla variabile T(t)
%%%%% dTi =DELAY(S×(1−GP),10)×(1−T)×0.1×AE  -----definizione di Stella
%
dTi = ((xd(6)*(exp(d*xd(5))))) * (1-x(7)) * 0.1 * AE;


f(7) = dTi - dTo ;

endfunction


t=[t_initial:delta_t:t_final];
res = ode45d (@fun, t, [init_l;init_f;init_m;init_u;init_lsd;init_s;init_t], [5*delta_t;10*delta_t],ones(7,10));

plot (res.x, res.y);
xlabel('Tempo (anni)')
ylabel('Variabili')
title('Patumod')
legend ("L", "F", "M", "U", "LSD", "S", "T")
axis([0,200, 0, 1])
replot
print('figureFRENK.ps', '-deps');

---

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

On 1 Nov 2009, at 20:08, franco basaglia wrote:

>
> I must be missing something as in your file I see no time  
> derivatives...
> sorry for my oversight. I attach again pdf with right system.

What exactly does "delta t" mean in your file?

> function f = fun (t, y, yd)
> f(1) =-y(1) + yd(2);
> f(2) =-y(2) ;
> endfunction
>
> ode23d (@fun, [0 5], [1;1], .1, [1;1])
>
> notice that in the example above only one equation includes a delay:
> is this what you mean by system of ODEs and DDEs?
>
> Yes,sure.
> To be more precise I have a system like this:
> d y1(t)/ dt = -y1(t)
> d y2(t)/ dt = -y2(t) + y1(t-5)
> d y3(t)/dt  = -y3(t) + y2(t-10)*y1(t-10)
>
> that I try to solve in this way:
> function f = fun (t, y, yd)
> f(1) =-y(1)
> f(2) =-y(2) + yd(1)
> f(3) =-y(3) + yd(2)*yd(1)
> endfunction
> t = [0:.5:20]
> res = ode45d (@fun, t, [1;1;1], [5;10], ones (3,10))

I think this should rather be something like:

function f = fun (t, y, yd)
f(1) =-y(1);                   %% y1' = -y1(t)
f(2) =-y(2) + yd(1,1);         %% y2' = -y2(t) + y1(t-lags(1))
f(3) =-y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) + y2(t-lags(2))*y1(t-
lags(2))
endfunction

T = [0,20]                                            %% only initlial  
and final time
                                                       %% need to be  
provided

res = ode45d (@fun, t, [1;1;1], [5, 10], ones (3,2)); %% number of  
columns in HIST should
                                                       %% be the same  
as the number of columns
                                                       %% in LAGS

I beleive Thomas Treichl, the author of odepkg is listening on the  
list so maybe he can confirm whether this is correct?

> It's run. But I don't think this is the right way.I don't know if  
> with this notation it applies a delay of 5 time units to f(2) and of  
> 10 time units to f(3).
what do you mean by "time units"?

> Note that in f(3) yd(1) comes again but with a different delay.
> Infact for my script, that I attach here, is the same. It's runs but  
> with different results compare to Stella.
>
> Any Idea?

 From your pdf file and your script I get the feeling that you are  
confusing the
time parameter in the continuos problem and the number of time-steps  
in the discretized
system...

> thanks a lot
>
> f.t.
c.

P.S. please keep the list in CC as someone else may want to help with  
your questions
or benefit from the answers.



[scriptode45.m]

---

%%%%valori iniziali delle coperture
%copertura iniziale di prato semplice(tra 0 e 1)
global init_l = 0.148;
%copertura iniziale di prato incolto(tra 0 e 1)
global init_f = 0.147;
%copertura iniziale di prato complesso(tra 0 e 1)
global init_m = 0.303;
%copertura iniziale di arbusti(tra 0 e 1)
global init_s = 0.123;
%copertura iniziale di alberi(tra 0 e 1)
global init_t = 0.402;
%copertura iniziale di sottobosco(tra 0 e 1)
global init_u = (1 - init_l - init_f - init_m);

%%%%variabili di controllo
%altitudine (m)
global A = 600;
%Carico Animale Globale (UBA giorni/ha/anno)
global  GSD =103;

%%%%parametri
%Valore Pastorale medio incolto
global PVF = 10;
%Valore Pastorale medio pascolo semplice
global PVL = 20;
%Valore Pastorale medio pascolo complesso
global PVM = 40;
%Valore Pastorale medio sottobosco
global PVU = 5;
%Valore Pastorale globale
global GPV = 21;

%%%%equazioni
%valore pastorale locale
global LPV = PVF*init_f+PVL*init_l+PVM*init_m+PVU*init_u;
%carico animale locale(UBA giorno/ha/anno)
global init_lsd = GSD*LPV/GPV ;
%parametri
global d = -0.013305;
global b = -0.005;
global c = 1081.093;
%effetti dell'altitudine(tra 0 e 1)
global AE = exp(b*(A-c))/(exp(b*(A-c))+1);
%capacità portante globale(UBA giorno/ha/anno)
global GCC = 115000*GPV/18/A;
%tasso di utilizzazione globale(tra 0 e 1)
global GU =GSD/GCC;


%%%%%%%% parametri temporali di simulazione
%t inizio simulazione
global t_initial = 0;
%t fine simulazione
global t_final = 200;
       
global monitor_points = 400;
global delta_t = (t_final - t_initial) / monitor_points;


function f = fun (t,x,xd)

global d
global GSD
global PVF
global PVL
global PVM
global PVU
global GPV
global b
global A
global c
global AE
global SCr
global SCs
global SCp
global TCr
global TCs
global TCp
global GCC
global GU
global init_t
global t_initial
global t_final
global monitor_points
global delta_t

%x1
global init_l  
%x2
global init_f
%x3
global init_m
%x4
global init_u
%x5
global LPV
%x6
global init_lsd
%x7
global init_s

%%%%%% L(t) = L(t−dt)+(−LtoM−LtoF)×dt Prato Semplice
%%%%%% L(t)= -((0.06*L*GP)-(0,05*M*(1-GP))) - ((0,5*L*(1-GP)*(O.3+S)*AE)-(0.1*F*GP))
f(1) = -((0.06*x(1)*(1-exp(d*x(5))))-(0.05*x(3)*(exp(d*x(5)))))-( (0.5*x(1)*(exp(d*x(5)))*(0.3+x(6)))*AE-(0.1*x(2)*(1-exp(d*x(5)))));

%%%%%% F(t) = F(t−dt)+(LtoF+MtoF−FtoU)×dt  Incolto
%%%%%% F(t)= ((0,5*L*(1-GP)*(O.3+S)*AE)-(0.1*F*GP)) + (O.8*M*(1-GP)*(O.3+S)-(0.4*F*(GP)^2) - (F*(T-U))
f(2) =(((0.5*x(1)*(exp(d*x(5)))*(0.3+x(6)))*AE-(0.1*x(2)*(1-exp(d*x(5))))) + (0.8*x(3)*(exp(d*x(5)))*(0.3+x(6))-(0.4*x(2)*((1-exp(d*x(5))^2)))) - (x(2)*(x(6)-x(4))));

%%%%%% M(t) = M(t−dt)+(LtoM−MtoF)×dt Prato complesso
%%%%%% M(t)=((0.06*L*GP)-(0,05*M*(1-GP))) - (O.8*M*(1-GP)*(O.3+S)-(0.4*F*(GP)^2)
f(3) = ((0.06*x(1)*(1-exp(d*x(5))))-(0.05*x(3)*(+exp(d*x(5))))) - (0.8*x(3)*(exp(d*x(5)))*(0.3+x(6))-(0.4*x(2)*((1-exp(d*x(5))^2))));

%%%%%% U(t) = U(t−dt)+(FtoU)×dt
%%%%%% U(t) = (F*(T-U))
f(4) = (x(2)*(x(6)-x(4)));

%%%%% if per calcolare LSD
LCC =LPV*115000/18/A;
LU = x(5)/LCC;
if (LU < 1)
  dLSD = x(5) * (GU-LU);
else
  dLSD = x(5) *(1-LU);
endif

%%%%%% LSD(t) =LSD(t−dt)+(dLSD)×dt capacità di carico locale
f(5) = dLSD ;



%%%%% S(t)=S(t−dt)+(dSi−dSo)×dt        Arbusti

dSo = (0.1 + 1 - exp(d*x(5))) * (0.1 + x(7)) * x(6) * 0.2;
%
%%%%% funzione DELAY sulla variabile S(t)
%%%%% dSi = 0.2×DELAY((F+0.5×U)×(1−GP),5)×(1−S)×(1−T)×AE -----definizione di Stella
%
dSi = 0.2*((xd(2)+0.5*xd(4))*(exp(d*xd(5))))*(1-x(6))*(1-x(7))*AE;


f(6) = dSi - dSo;



%%%%% T(t) = T(t−dt)+(dTi−dTo)×dt   Alberi

dTo = x(7)*0.005;
%
%%%%% funzione DELAY sulla variabile T(t)
%%%%% dTi =DELAY(S×(1−GP),10)×(1−T)×0.1×AE  -----definizione di Stella
%
dTi = ((xd(6)*(exp(d*xd(5))))) * (1-x(7)) * 0.1 * AE;


f(7) = dTi - dTo ;

endfunction


t=[t_initial:delta_t:t_final];
res = ode45d (@fun, t, [init_l;init_f;init_m;init_u;init_lsd;init_s;init_t], [5*delta_t;10*delta_t],ones(7,10));

plot (res.x, res.y);
xlabel('Tempo (anni)')
ylabel('Variabili')
title('Patumod')
legend ("L", "F", "M", "U", "LSD", "S", "T")
axis([0,200, 0, 1])
replot
print('figureFRENK.ps', '-deps');

---

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Thomas Treichl]

Carlo de Falco schrieb:

> On 1 Nov 2009, at 20:08, franco basaglia wrote:
>> Yes,sure.
>> To be more precise I have a system like this:
>> d y1(t)/ dt = -y1(t)
>> d y2(t)/ dt = -y2(t) + y1(t-5)
>> d y3(t)/dt  = -y3(t) + y2(t-10)*y1(t-10)
>>
>> that I try to solve in this way:
>> function f = fun (t, y, yd)
>> f(1) =-y(1)
>> f(2) =-y(2) + yd(1)
>> f(3) =-y(3) + yd(2)*yd(1)
>> endfunction
>> t = [0:.5:20]
>> res = ode45d (@fun, t, [1;1;1], [5;10], ones (3,10))
>
> I think this should rather be something like:
>
> function f = fun (t, y, yd)
> f(1) =-y(1);                   %% y1' = -y1(t)
> f(2) =-y(2) + yd(1,1);         %% y2' = -y2(t) + y1(t-lags(1))
> f(3) =-y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) +
> y2(t-lags(2))*y1(t-lags(2))
> endfunction
>
> T = [0,20]                                            %% only initlial
> and final time
>                                                       %% need to be
> provided
>
> res = ode45d (@fun, t, [1;1;1], [5, 10], ones (3,2)); %% number of
> columns in HIST should
>                                                       %% be the same as
> the number of columns
>                                                       %% in LAGS
> I beleive Thomas Treichl, the author of odepkg is listening on the list
> so maybe he can confirm whether this is correct?

Because there is an "ode" in the subject of the email, yes, but I also must say
that I currently don't have the time to dig into this problem and can only
provide a short answer: if the above equations form the DDE problem then Carlo's
implementation does make more sense to me.

Best regards

   Thomas
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[franco basaglia]

What exactly does "delta t" mean in your file?

what do you mean by "time units"?


From your pdf file and your script I get the feeling that you are confusing the
time parameter in the continuos problem and the number of time-steps in the discretized
system...

Yes.You're right. I have confused "delta t" (time parameter in the contiunuos problem) and time-steps in the discretizet system, that are my "time units".
These time units are 5 and 10 in my DISCRETE function "delay".
Are you meaning this or that I have to modify my system in a discrete time system?
I'm confuse about this question because, for example, Stella runtime equation is something like this:
....
L(t) = L(t−dt)+(−LtoM−LtoF)×dt
T(t) = T(t−dt)+(DELAY(S×(1−GP),10)×(1−T))×dt
and  Euler and 4th-order Runge-Kutta methods are used
with time step dt  set to 0.5 year


Anyway, I correct my script following your suggestion, but it doesn't work. And the different values from Stella concern, particularly, the 2 state variables with delay.
I have also this warning but I don't think that the problem is here:
warning: Option "RelTol" not set, new value 0.000001 is used
warning: Option "AbsTol" not set, new value 0.000001 is used
warning: Option "InitialStep" not set, new value 30.000000 is used
warning: Option "MaxStep" not set, new value 30.000000 is used


I attach my new script.



thanks

best regrads
f.b.

Il giorno 02 novembre 2009 00.00, Carlo de Falco <carlo.defalco@...> ha scritto:

On 1 Nov 2009, at 20:08, franco basaglia wrote:

I must be missing something as in your file I see no time derivatives...
sorry for my oversight. I attach again pdf with right system.

What exactly does "delta t" mean in your file?

function f = fun (t, y, yd)
f(1) =-y(1) + yd(2);
f(2) =-y(2) ;
endfunction

ode23d (@fun, [0 5], [1;1], .1, [1;1])

notice that in the example above only one equation includes a delay:
is this what you mean by system of ODEs and DDEs?

Yes,sure.
To be more precise I have a system like this:
d y1(t)/ dt = -y1(t)
d y2(t)/ dt = -y2(t) + y1(t-5)
d y3(t)/dt  = -y3(t) + y2(t-10)*y1(t-10)

that I try to solve in this way:
function f = fun (t, y, yd)
f(1) =-y(1)
f(2) =-y(2) + yd(1)
f(3) =-y(3) + yd(2)*yd(1)
endfunction
t = [0:.5:20]
res = ode45d (@fun, t, [1;1;1], [5;10], ones (3,10))

I think this should rather be something like:

function f = fun (t, y, yd)

f(1) =-y(1);                   %% y1' = -y1(t)
f(2) =-y(2) + yd(1,1);         %% y2' = -y2(t) + y1(t-lags(1))
f(3) =-y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) + y2(t-lags(2))*y1(t-lags(2))
endfunction

T = [0,20]                                            %% only initlial and final time
                                                     %% need to be provided

res = ode45d (@fun, t, [1;1;1], [5, 10], ones (3,2)); %% number of columns in HIST should
                                                     %% be the same as the number of columns
                                                     %% in LAGS

I beleive Thomas Treichl, the author of odepkg is listening on the list so maybe he can confirm whether this is correct?

It's run. But I don't think this is the right way.I don't know if with this notation it applies a delay of 5 time units to f(2) and of 10 time units to f(3).

what do you mean by "time units"?

Note that in f(3) yd(1) comes again but with a different delay.
Infact for my script, that I attach here, is the same. It's runs but with different results compare to Stella.

Any Idea?

>From your pdf file and your script I get the feeling that you are confusing the
time parameter in the continuos problem and the number of time-steps in the discretized
system...

thanks a lot

f.t.

c.

P.S. please keep the list in CC as someone else may want to help with your questions
or benefit from the answers.

[ode45delay.m]

---

%%%%valori iniziali delle coperture
%copertura iniziale di prato semplice(tra 0 e 1)
global init_l = 0.148;
%copertura iniziale di prato incolto(tra 0 e 1)
global init_f = 0.147;
%copertura iniziale di prato complesso(tra 0 e 1)
global init_m = 0.303;
%copertura iniziale di arbusti(tra 0 e 1)
global init_s = 0.123;
%copertura iniziale di alberi(tra 0 e 1)
global init_t = 0.402;
%copertura iniziale di sottobosco(tra 0 e 1)
global init_u = (1 - init_l - init_f - init_m);

%%%%variabili di controllo
%altitudine (m)
global A = 600;
%Carico Animale Globale (UBA giorni/ha/anno)
global  GSD =0;

%%%%parametri
%Valore Pastorale medio incolto
global PVF = 10;
%Valore Pastorale medio pascolo semplice
global PVL = 20;
%Valore Pastorale medio pascolo complesso
global PVM = 40;
%Valore Pastorale medio sottobosco
global PVU = 5;
%Valore Pastorale globale
global GPV = 21;

%%%%equazioni
%valore pastorale locale
global LPV = PVF*init_f+PVL*init_l+PVM*init_m+PVU*init_u;
%carico animale locale(UBA giorno/ha/anno)
global init_lsd = GSD*LPV/GPV ;
%parametri
global d = -0.013305;
global b = -0.005;
global c = 1081.093;
%effetti dell'altitudine(tra 0 e 1)
global AE = exp(b*(A-c))/(exp(b*(A-c))+1);
%capacità portante globale(UBA giorno/ha/anno)
global GCC = 115000*GPV/18/A;
%tasso di utilizzazione globale(tra 0 e 1)
global GU =GSD/GCC;


%%%%%%%% parametri temporali di simulazione
%t inizio simulazione
%global t_initial = 0;
%t fine simulazione
%global t_final = 300;
       
%global monitor_points = 200;
%global delta_t = (t_final - t_initial) / monitor_points;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = fun (t,x,xd)

global d
global GSD
global PVF
global PVL
global PVM
global PVU
global GPV
global b
global A
global c
global AE
global SCr
global SCs
global SCp
global TCr
global TCs
global TCp
global GCC
global GU
global init_t
global t_initial
global t_final
global monitor_points
global delta_t

%x1
global init_l  
%x2
global init_f
%x3
global init_m
%x4
global init_u
%x5
global LPV
%x6
global init_lsd
%x7
global init_s

%%%%%% L(t) = L(t−dt)+(−LtoM−LtoF)×dt Prato Semplice
%%%%%% L(t)= -((0.06*L*GP)-(0,05*M*(1-GP))) - ((0,5*L*(1-GP)*(O.3+S)*AE)-(0.1*F*GP))
f(1) = -((0.06*x(1)*(1-exp(d*x(5))))-(0.05*x(3)*(exp(d*x(5)))))-( (0.5*x(1)*(exp(d*x(5)))*(0.3+x(6)))*AE-(0.1*x(2)*(1-exp(d*x(5)))));

%%%%%% F(t) = F(t−dt)+(LtoF+MtoF−FtoU)×dt  Incolto
%%%%%% F(t)= ((0,5*L*(1-GP)*(O.3+S)*AE)-(0.1*F*GP)) + (O.8*M*(1-GP)*(O.3+S)-(0.4*F*(GP)^2) - (F*(T-U))
f(2) =(((0.5*x(1)*(exp(d*x(5)))*(0.3+x(6)))*AE-(0.1*x(2)*(1-exp(d*x(5))))) + (0.8*x(3)*(exp(d*x(5)))*(0.3+x(6))-(0.4*x(2)*((1-exp(d*x(5))^2)))) - (x(2)*(x(6)-x(4))));

%%%%%% M(t) = M(t−dt)+(LtoM−MtoF)×dt Prato complesso
%%%%%% M(t)=((0.06*L*GP)-(0,05*M*(1-GP))) - (O.8*M*(1-GP)*(O.3+S)-(0.4*F*(GP)^2)
f(3) = ((0.06*x(1)*(1-exp(d*x(5))))-(0.05*x(3)*(+exp(d*x(5))))) - (0.8*x(3)*(exp(d*x(5)))*(0.3+x(6))-(0.4*x(2)*((1-exp(d*x(5))^2))));

%%%%%% U(t) = U(t−dt)+(FtoU)×dt
%%%%%% U(t) = (F*(T-U))
f(4) = (x(2)*(x(6)-x(4)));

%%%%% if per calcolare LSD
LCC =LPV*115000/18/A;
LU = x(5)/LCC;
if (LU < 1)
  dLSD = x(5) * (GU-LU);
else
  dLSD = x(5) *(1-LU);
endif

%%%%%% LSD(t) =LSD(t−dt)+(dLSD)×dt capacità di carico locale
f(5) = dLSD ;



%%%%% S(t)=S(t−dt)+(dSi−dSo)×dt        Arbusti

dSo = (0.1 + 1 - exp(d*x(5))) * (0.1 + x(7)) * x(6) * 0.2;
%
%%%%% funzione DELAY sulla variabile S(t)
%%%%% dSi = 0.2×DELAY((F+0.5×U)×(1−GP),5)×(1−S)×(1−T)×AE -----definizione di Stella
%
dSi = 0.2*((xd(2,1)+0.5*xd(4,1))*(exp(d*xd(5,1))))*(1-x(6))*(1-x(7))*AE;


f(6) = dSi - dSo;



%%%%% T(t) = T(t−dt)+(dTi−dTo)×dt   Alberi

dTo = x(7)*0.005;
%
%%%%% funzione DELAY sulla variabile T(t)
%%%%% dTi =DELAY(S×(1−GP),10)×(1−T)×0.1×AE  -----definizione di Stella
%
dTi = ((xd(6,2)*(exp(d*xd(5,2))))) * (1-x(7)) * 0.1 * AE;


f(7) = dTi - dTo ;

endfunction


t=[0;300];


res = ode45d (@fun, t, [init_l;init_f;init_m;init_u;init_lsd;init_s;init_t], [5,10],ones(7,2));


plot (res.x, res.y);
xlabel('Tempo (anni)')
ylabel('Variabili')
title('Patumod')
legend ("L", "F", "M", "U", "LSD", "S", "T")
axis([0,200, 0, 1])
replot
print('figuregsd0a600.ps', '-deps');

---

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

On 3 Nov 2009, at 14:58, franco basaglia wrote:

>
>> What exactly does "delta t" mean in your file?
>>
>> what do you mean by "time units"?
>>
>>
>> From your pdf file and your script I get the feeling that you are  
>> confusing the
>> time parameter in the continuos problem and the number of time-
>> steps in the discretized
>> system...
>
> Yes.You're right. I have confused "delta t" (time parameter in the  
> contiunuos problem) and time-steps in the discretizet system, that  
> are my "time units".
> These time units are 5 and 10 in my DISCRETE function "delay".
> Are you meaning this or that I have to modify my system in a  
> discrete time system?
> I'm confuse about this question because, for example, Stella runtime  
> equation is something like this:
> ....
> L(t) = L(t−dt)+(−LtoM−LtoF)×dt
> T(t) = T(t−dt)+(DELAY(S×(1−GP),10)×(1−T))×dt
> and  Euler and 4th-order Runge-Kutta methods are used
> with time step dt  set to 0.5 year

the functions in odepkg use adaptive time-stepping so

1) you cannot know in advance what the value of dt will be

and

2) dt it will vary during the simulation

therefore specifying the delay in terms of dt makes no sense,
if your time is measured in years your delay must also be expressed in  
years.
In particular the following command in your script:

res = ode45d (@fun, t, init, [5,10],ones(7,2));

means that at any given time t, "fun" will be called as

fun (t, x, xd)

where xd(:,1) is the value of the state vetor at the time point t-5  
years
and xd(:,2) is the value of the state vetor at the time point t-10 years

is that what you want?

> Anyway, I correct my script following your suggestion, but it  
> doesn't work. And the different values from Stella concern,  
> particularly, the 2 state variables with delay.
> I have also this warning but I don't think that the problem is here:
> warning: Option "RelTol" not set, new value 0.000001 is used
> warning: Option "AbsTol" not set, new value 0.000001 is used
> warning: Option "InitialStep" not set, new value 30.000000 is used
> warning: Option "MaxStep" not set, new value 30.000000 is used

These are just warnings, you can safely ignore them at this stage, or  
specify the missing options if you want to get rid of them...

> I attach my new script.

Your script runs without errors on my system, but I have no idea about  
what the results should look like so I cannot comment on their  
correctness.

> thanks
>
> best regrads
> f.b.
c.
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[franco basaglia]

is that what you want?

Thanks Carlo. You are very helpful and I'm learning about Octave


Your script runs without errors on my system, but I have no idea about what the results should look like so I cannot comment on their correctness.

I make some editing on equations in my script(I attach it again).
I also try to run my system without delay, paradoxically it works better than "Delay Script".
I attach result from Stella( there isn't U variable), from Octave with Delay and From Octave without Delay.
Note that Variables S and T are influenced directly from Delay
As you can see the difference is on transitory.

Thanks for all.


Best regards
f.b.

2009/11/3 Carlo de Falco <carlo.defalco@...>

On 3 Nov 2009, at 14:58, franco basaglia wrote:

What exactly does "delta t" mean in your file?

what do you mean by "time units"?


>From your pdf file and your script I get the feeling that you are confusing the
time parameter in the continuos problem and the number of time-steps in the discretized
system...

Yes.You're right. I have confused "delta t" (time parameter in the contiunuos problem) and time-steps in the discretizet system, that are my "time units".
These time units are 5 and 10 in my DISCRETE function "delay".
Are you meaning this or that I have to modify my system in a discrete time system?
I'm confuse about this question because, for example, Stella runtime equation is something like this:
....
L(t) = L(t−dt)+(−LtoM−LtoF)×dt
T(t) = T(t−dt)+(DELAY(S×(1−GP),10)×(1−T))×dt
and  Euler and 4th-order Runge-Kutta methods are used
with time step dt  set to 0.5 year

the functions in odepkg use adaptive time-stepping so

1) you cannot know in advance what the value of dt will be

and

2) dt it will vary during the simulation

therefore specifying the delay in terms of dt makes no sense,
if your time is measured in years your delay must also be expressed in years.
In particular the following command in your script:

res = ode45d (@fun, t, init, [5,10],ones(7,2));

means that at any given time t, "fun" will be called as

fun (t, x, xd)

where xd(:,1) is the value of the state vetor at the time point t-5 years
and xd(:,2) is the value of the state vetor at the time point t-10 years

is that what you want?

Anyway, I correct my script following your suggestion, but it doesn't work. And the different values from Stella concern, particularly, the 2 state variables with delay.
I have also this warning but I don't think that the problem is here:
warning: Option "RelTol" not set, new value 0.000001 is used
warning: Option "AbsTol" not set, new value 0.000001 is used
warning: Option "InitialStep" not set, new value 30.000000 is used
warning: Option "MaxStep" not set, new value 30.000000 is used

These are just warnings, you can safely ignore them at this stage, or specify the missing options if you want to get rid of them...

I attach my new script.

Your script runs without errors on my system, but I have no idea about what the results should look like so I cannot comment on their correctness.

thanks

best regrads
f.b.

c.

[patumode45d.m]

---

%%%-------valori iniziali delle coperture
%--copertura iniziale di prato semplice(tra 0 e 1)
global init_L = 0.148;
%--copertura iniziale di prato incolto(tra 0 e 1)
global init_F = 0.147;
%--copertura iniziale di prato complesso(tra 0 e 1)
global init_M = 0.303;
%--copertura iniziale di arbusti(tra 0 e 1)
global init_S = 0.123;
%--copertura iniziale di alberi(tra 0 e 1)
global init_T = 0.402;
%--copertura iniziale di sottobosco(tra 0 e 1)
global init_U = (1 - init_L - init_F - init_M);


%%%-------variabili di controllo
%--altitudine (m)
global A = 600;
%--Carico Animale Globale (UBA giorni/ha/anno)
global  GSD =0;


%%%%-------parametri
%--Valore Pastorale medio incolto
global PVF = 10;
%--Valore Pastorale medio pascolo semplice
global PVL = 20;
%--Valore Pastorale medio pascolo complesso
global PVM = 40;
%--Valore Pastorale medio sottobosco
global PVU = 5;
%--Valore Pastorale globale
global GPV = 21;


%%%%-------relazioni
%--valore pastorale locale
global init_LPV = 1*PVF*init_F+PVL*init_L+PVM*init_M+1*PVU*init_U;
%--carico animale locale (UBA giorno/ha/anno)
global init_LSD = GSD*init_LPV/GPV ;
%--parametri
global d = -0.013305;
global b = -0.005;
global c = 1081.093;
%--effetti dell'altitudine (tra 0 e 1)
global AE = exp(b*(A-c))/(exp(b*(A-c))+1);
%--capacità portante globale(UBA giorno/ha/anno)
global GCC = 115000*GPV/18/A;
%--tasso di utilizzazione globale(tra 0 e 1)
global GU = GSD/GCC;
%--capacità portante locale (UBA giorno/ha/anno)
global init_LCC = init_LPV*115000/18/A;
%--tasso di utilizzazione locale (tra 0 e 1)
global init_LU = init_LSD/init_LCC

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = fun (t,x,xd)

global GSD
global A
global d
global b
global c
global PVF
global PVL
global PVM
global PVU
global GPV
global AE
%global SCr
%global SCs
%global SCp
%global TCr
%global TCs
%global TCp
global GCC
global GU
global LPV
global LCC
global LU
global init_LU

%x1
global init_L
%x2
global init_F
%x3
global init_M
%x4
global init_U
%x5
global init_LSD
%x6
global init_S
%x7
global init_T





%%%%%% L(t) = L(t−dt)+(−LtoM−LtoF)×dt Prato Semplice
%%%%%% L(t)= -((0.06*L*GP)-(0.05*M*(1-GP))) - ((0.5*L*(1-GP)*(O.3+S)*AE)-(0.1*F*GP))
f(1) =  - ((0.06*x(1)*(1-exp(d*x(5))))-(0.05*x(3)*exp(d*x(5)))) - (((0.5*x(1)*(exp(d*x(5)))*(0.3+x(6))*AE))-(0.1*x(2)*(1-exp(d*x(5)))));

%%%%%% F(t) = F(t−dt)+(LtoF+MtoF−FtoU)×dt  Incolto
%%%%%% F(t)= ((0,5*L*(1-GP)*(O.3+S)*AE)-(0.1*F*GP)) + (O.8*M*(1-GP)*(O.3+S)-(0.4*F*(GP)^2) - (F*(T-U))
f(2)= (((0.5*x(1)*(exp(d*x(5)))*(0.3+x(6))*AE))-(0.1*x(2)*(1-exp(d*x(5))))) + ((0.8*x(3)*exp(d*x(5)))*(0.3+x(6))-(0.4*x(2)*(1-exp(d*x(5)))^2)) - (x(2)*(x(7)-x(4)));


%%%%%% M(t) = M(t−dt)+(LtoM−MtoF)×dt Prato complesso
%%%%%% M(t)=((0.06*L*GP)-(0,05*M*(1-GP))) - (O.8*M*(1-GP)*(O.3+S)-(0.4*F*(GP)^2)
f(3)= ((0.06*x(1)*(1-exp(d*x(5))))-(0.05*x(3)*exp(d*x(5)))) - ((0.8*x(3)*exp(d*x(5)))*(0.3+x(6))-(0.4*x(2)*(1-exp(d*x(5)))^2));


%%%%%% U(t) = U(t−dt)+(FtoU)×dt
%%%%%% U(t) = (F*(T-U))
f(4)= (x(2)*(x(7)-x(4)));



%%%%% if per calcolare LSD


if (t=0)
        if (init_LU < 1)
          dLSD = x(5) * (GU-LU);
        else
          dLSD = x(5) *(1-LU);
        endif
else
        if (LU < 1)
          dLSD = x(5) * (GU-LU);
        else
          dLSD = x(5) *(1-LU);
        endif
endif

%%%%%% LSD(t) =LSD(t−dt)+(dLSD)×dt capacità di carico locale
f(5) = dLSD;


%%%%% S(t)=S(t−dt)+(dSi−dSo)×dt        Arbusti
%%%%% S(t) = (0.2×DELAY((F+0.5×U)×(1−GP),5)×(1−S)×(1−T)×AE)-((0.1+GP)*(0.1+T)*S*0.2)---definizione del delay secondo Stella
f(6) =(0.2*((xd(2,1)+0.5*xd(4,1))*(exp(d*xd(5,1))))*(1-x(6))*(1-x(7))*AE) - ((0.1 +(1-exp(d*x(5))))*(0.1+x(7))*x(6)* 0.2);


%%%%% T(t) = T(t−dt)+(dTi−dTo)×dt   Alberi
%%%%% T(t) = (DELAY(S×(1−GP),10)×(1−T)×0.1×AE) - T*0.005 ---definizione del delay secondo Stella
f(7) = ((xd(6,2)*(exp(d*xd(5,2))))*(1-x(7))*0.1*AE) - (x(7)*0.005);


LPV = 1*PVF*x(2)+PVL*x(1)+PVM*x(3)+1*PVU*x(4)

LCC = LPV*115000/18/A

LU = x(5)/LCC


endfunction


t=[0:200];


res = ode45d (@fun, t, [init_L;init_F;init_M;init_U;init_LSD;init_S;init_T], [5,10],ones(7,2));


plot (res.x, res.y);
xlabel('Tempo (anni)')
ylabel('Variabili')
title('Patumod')
legend ("L", "F", "M", "U", "LSD", "S", "T")
axis([0,100, 0, 1])
replot
print('OctDelaygsd0a600.ps', '-deps');

---

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

Thomas,

On 2 Nov 2009, at 20:26, Thomas Treichl wrote:

> Because there is an "ode" in the subject of the email, yes, but I  
> also must say that I currently don't have the time to dig into this  
> problem and can only provide a short answer: if the above equations  
> form the DDE problem then Carlo's implementation does make more  
> sense to me.

Thanks for your answer, I was not 100% sure I was doing things  
correctly as
it seems to me that the case of mutiple delays is not fully clear from  
the documentation (I looked at the implementation to find out the  
meaning of LAGS).
So if you believe my interpretation was correct, would you mind if I  
commit the patch below to improve the documentation of the DDE  
functions?

> Best regards
>  Thomas

c.

Index: ode23d.m
===================================================================
--- ode23d.m (revision 6430)
+++ ode23d.m (working copy)
@@ -23,6 +23,14 @@
 %#
 %# If this function is called with no return argument then plot the solution over time in a figure window while solving the set of DDEs that are defined in a function and specified by the function handle @var{@@fun}. The second input argument @var{slot} is a double vector that defines the time slot, @var{init} is a double vector that defines the initial values of the states, @var{lags} is a double vector that describes the lags of time, @var{hist} is a double matrix and describes the history of the DDEs, @var{opt} can optionally be a structure array that keeps the options created with the command @command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be other input arguments of any type that have to be passed to the function defined by @var{@@fun}.
 %#
+%# In other words, this function will solve a problem of the form
+%#
+%# @example
+%# dy/dt = fun (t, y(t), y(t-lags(1), y(t-lags(2), @dots{})))
+%# y(slot(1)) = init
+%# y(slot(1)-lags(1)) = hist(1), y(slot(1)-lags(2)) = hist(2), @dots{}
+%# @end example
+%#
 %# If this function is called with one return argument then return the solution @var{sol} of type structure array after solving the set of DDEs. The solution @var{sol} has the fields @var{x} of type double column vector for the steps chosen by the solver, @var{y} of type double column vector for the solutions at each time step of @var{x}, @var{solver} of type string for the solver name and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector that keep the informations of the event function if an event function handle is set in the option argument @var{opt}.
 %#
 %# If this function is called with more than one return argument then return the time stamps @var{t}, the solution values @var{y} and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector.
@@ -37,6 +45,21 @@
 %# vlag = interp1 (vsol.x, vsol.y, vsol.x - 2);
 %# plot (vsol.y, vlag); legend ("fcao (t,y,z)");
 %# @end example
+%# Solve the following problem with two delayed state variables
+%# @example
+%# d y1(t)/ dt = -y1(t)
+%# d y2(t)/ dt = -y2(t) + y1(t-5)
+%# d y3(t)/dt  = -y3(t) + y2(t-10)*y1(t-10)
+%# @end example
+%# @example
+%# function f = fun (t, y, yd)
+%# f(1) =-y(1);                   %% y1' = -y1(t)
+%# f(2) =-y(2) + yd(1,1);         %% y2' = -y2(t) + y1(t-lags(1))
+%# f(3) =-y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) + y2(t-lags(2))*y1(t-lags(2))
+%# endfunction
+%# T = [0,20]                                                                                            
+%# res = ode45d (@fun, t, [1;1;1], [5, 10], ones (3,2));                                                  
+%# @end example
 %# @end deftypefn
 %#
 %# @seealso{odepkg}

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Thomas Treichl]

Carlo de Falco schrieb:

> Thomas,
>
> On 2 Nov 2009, at 20:26, Thomas Treichl wrote:
>
>> Because there is an "ode" in the subject of the email, yes, but I also
>> must say that I currently don't have the time to dig into this problem
>> and can only provide a short answer: if the above equations form the
>> DDE problem then Carlo's implementation does make more sense to me.
>
> Thanks for your answer, I was not 100% sure I was doing things correctly as
> it seems to me that the case of mutiple delays is not fully clear from
> the documentation (I looked at the implementation to find out the
> meaning of LAGS).
> So if you believe my interpretation was correct, would you mind if I
> commit the patch below to improve the documentation of the DDE functions?

These changes would be an improvement.
Can you please make these changes also for ode45d, ode54d and ode78d?

Thanks, Thomas
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

On 5 Nov 2009, at 16:48, Thomas Treichl wrote:

> Carlo de Falco schrieb:
>> Thomas,
>> On 2 Nov 2009, at 20:26, Thomas Treichl wrote:
>>> Because there is an "ode" in the subject of the email, yes, but I  
>>> also must say that I currently don't have the time to dig into  
>>> this problem and can only provide a short answer: if the above  
>>> equations form the DDE problem then Carlo's implementation does  
>>> make more sense to me.
>> Thanks for your answer, I was not 100% sure I was doing things  
>> correctly as
>> it seems to me that the case of mutiple delays is not fully clear  
>> from the documentation (I looked at the implementation to find out  
>> the meaning of LAGS).
>> So if you believe my interpretation was correct, would you mind if  
>> I commit the patch below to improve the documentation of the DDE  
>> functions?
>
> These changes would be an improvement.
> Can you please make these changes also for ode45d, ode54d and ode78d?
>
> Thanks, Thomas

I checked-in the patch below,
c.

Index: ode45d.m
===================================================================
--- ode45d.m (revision 6438)
+++ ode45d.m (revision 6439)
@@ -23,20 +23,54 @@
  %#
  %# If this function is called with no return argument then plot the  
solution over time in a figure window while solving the set of DDEs  
that are defined in a function and specified by the function handle  
@var{@@fun}. The second input argument @var{slot} is a double vector  
that defines the time slot, @var{init} is a double vector that defines  
the initial values of the states, @var{lags} is a double vector that  
describes the lags of time, @var{hist} is a double matrix and  
describes the history of the DDEs, @var{opt} can optionally be a  
structure array that keeps the options created with the command  
@command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be  
other input arguments of any type that have to be passed to the  
function defined by @var{@@fun}.
  %#
+%# In other words, this function will solve a problem of the form
+%# @example
+%# dy/dt = fun (t, y(t), y(t-lags(1), y(t-lags(2), @dots{})))
+%# y(slot(1)) = init
+%# y(slot(1)-lags(1)) = hist(1), y(slot(1)-lags(2)) = hist(2), @dots{}
+%# @end example
+%#
  %# If this function is called with one return argument then return  
the solution @var{sol} of type structure array after solving the set  
of DDEs. The solution @var{sol} has the fields @var{x} of type double  
column vector for the steps chosen by the solver, @var{y} of type  
double column vector for the solutions at each time step of @var{x},  
@var{solver} of type string for the solver name and optionally the  
extended time stamp information @var{xe}, the extended solution  
information @var{ye} and the extended index information @var{ie} all  
of type double column vector that keep the informations of the event  
function if an event function handle is set in the option argument  
@var{opt}.
  %#
  %# If this function is called with more than one return argument  
then return the time stamps @var{t}, the solution values @var{y} and  
optionally the extended time stamp information @var{xe}, the extended  
solution information @var{ye} and the extended index information  
@var{ie} all of type double column vector.
  %#
-%# For example, solve an anonymous implementation of a chaotic behavior
+%# For example:
+%# @itemize @minus
+%# @item
+%# the following code solves an anonymous implementation of a chaotic  
behavior
+%#
  %# @example
  %# fcao = @@(vt, vy, vz) [2 * vz / (1 + vz^9.65) - vy];
  %#
-%# vopt = odeset ("NormControl", "on", "RelTol", 1e-4);
+%# vopt = odeset ("NormControl", "on", "RelTol", 1e-3);
  %# vsol = ode45d (fcao, [0, 100], 0.5, 2, 0.5, vopt);
  %#
  %# vlag = interp1 (vsol.x, vsol.y, vsol.x - 2);
  %# plot (vsol.y, vlag); legend ("fcao (t,y,z)");
  %# @end example
+%#
+%# @item
+%# to solve the following problem with two delayed state variables
+%#
+%# @example
+%# d y1(t)/ dt = -y1(t)
+%# d y2(t)/ dt = -y2(t) + y1(t-5)
+%# d y3(t)/dt  = -y3(t) + y2(t-10)*y1(t-10)
+%# @end example
+%#
+%# one might do the following
+%#
+%# @example
+%# function f = fun (t, y, yd)
+%# f(1) =-y(1);                   %% y1' = -y1(t)
+%# f(2) =-y(2) + yd(1,1);         %% y2' = -y2(t) + y1(t-lags(1))
+%# f(3) =-y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) + y2(t-
lags(2))*y1(t-lags(2))
+%# endfunction
+%# T = [0,20]
+%# res = ode45d (@fun, t, [1;1;1], [5, 10], ones (3,2));
+%# @end example
+%#
+%# @end itemize
  %# @end deftypefn
  %#
  %# @seealso{odepkg}
Index: ode54d.m
===================================================================
--- ode54d.m (revision 6438)
+++ ode54d.m (revision 6439)
@@ -23,11 +23,21 @@
  %#
  %# If this function is called with no return argument then plot the  
solution over time in a figure window while solving the set of DDEs  
that are defined in a function and specified by the function handle  
@var{@@fun}. The second input argument @var{slot} is a double vector  
that defines the time slot, @var{init} is a double vector that defines  
the initial values of the states, @var{lags} is a double vector that  
describes the lags of time, @var{hist} is a double matrix and  
describes the history of the DDEs, @var{opt} can optionally be a  
structure array that keeps the options created with the command  
@command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be  
other input arguments of any type that have to be passed to the  
function defined by @var{@@fun}.
  %#
+%# In other words, this function will solve a problem of the form
+%# @example
+%# dy/dt = fun (t, y(t), y(t-lags(1), y(t-lags(2), @dots{})))
+%# y(slot(1)) = init
+%# y(slot(1)-lags(1)) = hist(1), y(slot(1)-lags(2)) = hist(2), @dots{}
+%# @end example
+%#
  %# If this function is called with one return argument then return  
the solution @var{sol} of type structure array after solving the set  
of DDEs. The solution @var{sol} has the fields @var{x} of type double  
column vector for the steps chosen by the solver, @var{y} of type  
double column vector for the solutions at each time step of @var{x},  
@var{solver} of type string for the solver name and optionally the  
extended time stamp information @var{xe}, the extended solution  
information @var{ye} and the extended index information @var{ie} all  
of type double column vector that keep the informations of the event  
function if an event function handle is set in the option argument  
@var{opt}.
  %#
  %# If this function is called with more than one return argument  
then return the time stamps @var{t}, the solution values @var{y} and  
optionally the extended time stamp information @var{xe}, the extended  
solution information @var{ye} and the extended index information  
@var{ie} all of type double column vector.
  %#
-%# For example, solve an anonymous implementation of a chaotic behavior
+%# For example:
+%# @itemize @minus
+%# @item
+%# the following code solves an anonymous implementation of a chaotic  
behavior
  %# @example
  %# fcao = @@(vt, vy, vz) [2 * vz / (1 + vz^9.65) - vy];
  %#
@@ -37,6 +47,29 @@
  %# vlag = interp1 (vsol.x, vsol.y, vsol.x - 2);
  %# plot (vsol.y, vlag); legend ("fcao (t,y,z)");
  %# @end example
+%#
+%# @item
+%# to solve the following problem with two delayed state variables
+%#
+%# @example
+%# d y1(t)/ dt = -y1(t)
+%# d y2(t)/ dt = -y2(t) + y1(t-5)
+%# d y3(t)/dt  = -y3(t) + y2(t-10)*y1(t-10)
+%# @end example
+%#
+%# one might do the following
+%#
+%# @example
+%# function f = fun (t, y, yd)
+%# f(1) =-y(1);                   %% y1' = -y1(t)
+%# f(2) =-y(2) + yd(1,1);         %% y2' = -y2(t) + y1(t-lags(1))
+%# f(3) =-y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) + y2(t-
lags(2))*y1(t-lags(2))
+%# endfunction
+%# T = [0,20]
+%# res = ode54d (@fun, t, [1;1;1], [5, 10], ones (3,2));
+%# @end example
+%#
+%# @end itemize
  %# @end deftypefn
  %#
  %# @seealso{odepkg}
Index: ode78d.m
===================================================================
--- ode78d.m (revision 6438)
+++ ode78d.m (revision 6439)
@@ -23,11 +23,21 @@
  %#
  %# If this function is called with no return argument then plot the  
solution over time in a figure window while solving the set of DDEs  
that are defined in a function and specified by the function handle  
@var{@@fun}. The second input argument @var{slot} is a double vector  
that defines the time slot, @var{init} is a double vector that defines  
the initial values of the states, @var{lags} is a double vector that  
describes the lags of time, @var{hist} is a double matrix and  
describes the history of the DDEs, @var{opt} can optionally be a  
structure array that keeps the options created with the command  
@command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be  
other input arguments of any type that have to be passed to the  
function defined by @var{@@fun}.
  %#
+%# In other words, this function will solve a problem of the form
+%# @example
+%# dy/dt = fun (t, y(t), y(t-lags(1), y(t-lags(2), @dots{})))
+%# y(slot(1)) = init
+%# y(slot(1)-lags(1)) = hist(1), y(slot(1)-lags(2)) = hist(2), @dots{}
+%# @end example
+%#
  %# If this function is called with one return argument then return  
the solution @var{sol} of type structure array after solving the set  
of DDEs. The solution @var{sol} has the fields @var{x} of type double  
column vector for the steps chosen by the solver, @var{y} of type  
double column vector for the solutions at each time step of @var{x},  
@var{solver} of type string for the solver name and optionally the  
extended time stamp information @var{xe}, the extended solution  
information @var{ye} and the extended index information @var{ie} all  
of type double column vector that keep the informations of the event  
function if an event function handle is set in the option argument  
@var{opt}.
  %#
  %# If this function is called with more than one return argument  
then return the time stamps @var{t}, the solution values @var{y} and  
optionally the extended time stamp information @var{xe}, the extended  
solution information @var{ye} and the extended index information  
@var{ie} all of type double column vector.
  %#
-%# For example, solve an anonymous implementation of a chaotic behavior
+%# For example:
+%# @itemize @minus
+%# @item
+%# the following code solves an anonymous implementation of a chaotic  
behavior
  %# @example
  %# fcao = @@(vt, vy, vz) [2 * vz / (1 + vz^9.65) - vy];
  %#
@@ -37,6 +47,29 @@
  %# vlag = interp1 (vsol.x, vsol.y, vsol.x - 2);
  %# plot (vsol.y, vlag); legend ("fcao (t,y,z)");
  %# @end example
+%#
+%# @item
+%# to solve the following problem with two delayed state variables
+%#
+%# @example
+%# d y1(t)/ dt = -y1(t)
+%# d y2(t)/ dt = -y2(t) + y1(t-5)
+%# d y3(t)/dt  = -y3(t) + y2(t-10)*y1(t-10)
+%# @end example
+%#
+%# one might do the following
+%#
+%# @example
+%# function f = fun (t, y, yd)
+%# f(1) =-y(1);                   %% y1' = -y1(t)
+%# f(2) =-y(2) + yd(1,1);         %% y2' = -y2(t) + y1(t-lags(1))
+%# f(3) =-y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) + y2(t-
lags(2))*y1(t-lags(2))
+%# endfunction
+%# T = [0,20]
+%# res = ode78d (@fun, t, [1;1;1], [5, 10], ones (3,2));
+%# @end example
+%#
+%# @end itemize
  %# @end deftypefn
  %#
  %# @seealso{odepkg}
Index: ode23d.m
===================================================================
--- ode23d.m (revision 6438)
+++ ode23d.m (revision 6439)
@@ -23,11 +23,22 @@
  %#
  %# If this function is called with no return argument then plot the  
solution over time in a figure window while solving the set of DDEs  
that are defined in a function and specified by the function handle  
@var{@@fun}. The second input argument @var{slot} is a double vector  
that defines the time slot, @var{init} is a double vector that defines  
the initial values of the states, @var{lags} is a double vector that  
describes the lags of time, @var{hist} is a double matrix and  
describes the history of the DDEs, @var{opt} can optionally be a  
structure array that keeps the options created with the command  
@command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be  
other input arguments of any type that have to be passed to the  
function defined by @var{@@fun}.
  %#
+%# In other words, this function will solve a problem of the form
+%# @example
+%# dy/dt = fun (t, y(t), y(t-lags(1), y(t-lags(2), @dots{})))
+%# y(slot(1)) = init
+%# y(slot(1)-lags(1)) = hist(1), y(slot(1)-lags(2)) = hist(2), @dots{}
+%# @end example
+%#
  %# If this function is called with one return argument then return  
the solution @var{sol} of type structure array after solving the set  
of DDEs. The solution @var{sol} has the fields @var{x} of type double  
column vector for the steps chosen by the solver, @var{y} of type  
double column vector for the solutions at each time step of @var{x},  
@var{solver} of type string for the solver name and optionally the  
extended time stamp information @var{xe}, the extended solution  
information @var{ye} and the extended index information @var{ie} all  
of type double column vector that keep the informations of the event  
function if an event function handle is set in the option argument  
@var{opt}.
  %#
  %# If this function is called with more than one return argument  
then return the time stamps @var{t}, the solution values @var{y} and  
optionally the extended time stamp information @var{xe}, the extended  
solution information @var{ye} and the extended index information  
@var{ie} all of type double column vector.
  %#
-%# For example, solve an anonymous implementation of a chaotic behavior
+%# For example:
+%# @itemize @minus
+%# @item
+%# the following code solves an anonymous implementation of a chaotic  
behavior
+%#
  %# @example
  %# fcao = @@(vt, vy, vz) [2 * vz / (1 + vz^9.65) - vy];
  %#
@@ -37,6 +48,29 @@
  %# vlag = interp1 (vsol.x, vsol.y, vsol.x - 2);
  %# plot (vsol.y, vlag); legend ("fcao (t,y,z)");
  %# @end example
+%#
+%# @item
+%# to solve the following problem with two delayed state variables
+%#
+%# @example
+%# d y1(t)/ dt = -y1(t)
+%# d y2(t)/ dt = -y2(t) + y1(t-5)
+%# d y3(t)/dt  = -y3(t) + y2(t-10)*y1(t-10)
+%# @end example
+%#
+%# one might do the following
+%#
+%# @example
+%# function f = fun (t, y, yd)
+%# f(1) =-y(1);                   %% y1' = -y1(t)
+%# f(2) =-y(2) + yd(1,1);         %% y2' = -y2(t) + y1(t-lags(1))
+%# f(3) =-y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) + y2(t-
lags(2))*y1(t-lags(2))
+%# endfunction
+%# T = [0,20]
+%# res = ode23d (@fun, t, [1;1;1], [5, 10], ones (3,2));
+%# @end example
+%#
+%# @end itemize
  %# @end deftypefn
  %#
  %# @seealso{odepkg}

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[franco basaglia]

Thank you all for your work.

I have a last question.
I'd convert my ode system to a discrete-time system
with annual time step. In thsi way I can simplify my delay question.

 What is the octave solver to use? How can I plan the system in Octave?

The discrete system is something like that:

y1(t) = y1(t-1) -y1(t)
y2(t) = y2(t-1) -y2(t) + y1(t-5)
y3(t) = t3(t-1) -y3(t) + y2(t-10)*y1(t-10)



Thank you

best regards
f.b.

2009/11/5 Carlo de Falco <carlo.defalco@...>

On 5 Nov 2009, at 16:48, Thomas Treichl wrote:

Carlo de Falco schrieb:

Thomas,
On 2 Nov 2009, at 20:26, Thomas Treichl wrote:

Because there is an "ode" in the subject of the email, yes, but I also must say that I currently don't have the time to dig into this problem and can only provide a short answer: if the above equations form the DDE problem then Carlo's implementation does make more sense to me.

Thanks for your answer, I was not 100% sure I was doing things correctly as
it seems to me that the case of mutiple delays is not fully clear from the documentation (I looked at the implementation to find out the meaning of LAGS).
So if you believe my interpretation was correct, would you mind if I commit the patch below to improve the documentation of the DDE functions?

These changes would be an improvement.
Can you please make these changes also for ode45d, ode54d and ode78d?

Thanks, Thomas

I checked-in the patch below,
c.

Index: ode45d.m
===================================================================
--- ode45d.m    (revision 6438)
+++ ode45d.m    (revision 6439)
@@ -23,20 +23,54 @@
 %#
 %# If this function is called with no return argument then plot the solution over time in a figure window while solving the set of DDEs that are defined in a function and specified by the function handle @var{@@fun}. The second input argument @var{slot} is a double vector that defines the time slot, @var{init} is a double vector that defines the initial values of the states, @var{lags} is a double vector that describes the lags of time, @var{hist} is a double matrix and describes the history of the DDEs, @var{opt} can optionally be a structure array that keeps the options created with the command @command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be other input arguments of any type that have to be passed to the function defined by @var{@@fun}.
 %#
+%# In other words, this function will solve a problem of the form
+%# @example
+%# dy/dt = fun (t, y(t), y(t-lags(1), y(t-lags(2), @dots{})))
+%# y(slot(1)) = init
+%# y(slot(1)-lags(1)) = hist(1), y(slot(1)-lags(2)) = hist(2), @dots{}
+%# @end example
+%#
 %# If this function is called with one return argument then return the solution @var{sol} of type structure array after solving the set of DDEs. The solution @var{sol} has the fields @var{x} of type double column vector for the steps chosen by the solver, @var{y} of type double column vector for the solutions at each time step of @var{x}, @var{solver} of type string for the solver name and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector that keep the informations of the event function if an event function handle is set in the option argument @var{opt}.
 %#
 %# If this function is called with more than one return argument then return the time stamps @var{t}, the solution values @var{y} and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector.
 %#
-%# For example, solve an anonymous implementation of a chaotic behavior
+%# For example:
+%# @itemize @minus
+%# @item
+%# the following code solves an anonymous implementation of a chaotic behavior
+%#
 %# @example
 %# fcao = @@(vt, vy, vz) [2 * vz / (1 + vz^9.65) - vy];
 %#
-%# vopt = odeset ("NormControl", "on", "RelTol", 1e-4);
+%# vopt = odeset ("NormControl", "on", "RelTol", 1e-3);
 %# vsol = ode45d (fcao, [0, 100], 0.5, 2, 0.5, vopt);
 %#
 %# vlag = interp1 (vsol.x, vsol.y, vsol.x - 2);
 %# plot (vsol.y, vlag); legend ("fcao (t,y,z)");
 %# @end example
+%#
+%# @item
+%# to solve the following problem with two delayed state variables
+%#
+%# @example
+%# d y1(t)/ dt = -y1(t)

+%# d y2(t)/ dt = -y2(t) + y1(t-5)
+%# d y3(t)/dt  = -y3(t) + y2(t-10)*y1(t-10)

+%# @end example
+%#
+%# one might do the following
+%#
+%# @example
+%# function f = fun (t, y, yd)

+%# f(1) =-y(1);                   %% y1' = -y1(t)
+%# f(2) =-y(2) + yd(1,1);         %% y2' = -y2(t) + y1(t-lags(1))
+%# f(3) =-y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) + y2(t-lags(2))*y1(t-lags(2))
+%# endfunction
+%# T = [0,20]

+%# res = ode45d (@fun, t, [1;1;1], [5, 10], ones (3,2));
+%# @end example
+%#
+%# @end itemize
 %# @end deftypefn
 %#
 %# @seealso{odepkg}
Index: ode54d.m
===================================================================
--- ode54d.m    (revision 6438)
+++ ode54d.m    (revision 6439)
@@ -23,11 +23,21 @@
 %#
 %# If this function is called with no return argument then plot the solution over time in a figure window while solving the set of DDEs that are defined in a function and specified by the function handle @var{@@fun}. The second input argument @var{slot} is a double vector that defines the time slot, @var{init} is a double vector that defines the initial values of the states, @var{lags} is a double vector that describes the lags of time, @var{hist} is a double matrix and describes the history of the DDEs, @var{opt} can optionally be a structure array that keeps the options created with the command @command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be other input arguments of any type that have to be passed to the function defined by @var{@@fun}.
 %#
+%# In other words, this function will solve a problem of the form
+%# @example
+%# dy/dt = fun (t, y(t), y(t-lags(1), y(t-lags(2), @dots{})))
+%# y(slot(1)) = init
+%# y(slot(1)-lags(1)) = hist(1), y(slot(1)-lags(2)) = hist(2), @dots{}
+%# @end example
+%#
 %# If this function is called with one return argument then return the solution @var{sol} of type structure array after solving the set of DDEs. The solution @var{sol} has the fields @var{x} of type double column vector for the steps chosen by the solver, @var{y} of type double column vector for the solutions at each time step of @var{x}, @var{solver} of type string for the solver name and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector that keep the informations of the event function if an event function handle is set in the option argument @var{opt}.
 %#
 %# If this function is called with more than one return argument then return the time stamps @var{t}, the solution values @var{y} and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector.
 %#
-%# For example, solve an anonymous implementation of a chaotic behavior
+%# For example:
+%# @itemize @minus
+%# @item
+%# the following code solves an anonymous implementation of a chaotic behavior
 %# @example
 %# fcao = @@(vt, vy, vz) [2 * vz / (1 + vz^9.65) - vy];
 %#
@@ -37,6 +47,29 @@
 %# vlag = interp1 (vsol.x, vsol.y, vsol.x - 2);
 %# plot (vsol.y, vlag); legend ("fcao (t,y,z)");
 %# @end example
+%#
+%# @item
+%# to solve the following problem with two delayed state variables
+%#
+%# @example
+%# d y1(t)/ dt = -y1(t)

+%# d y2(t)/ dt = -y2(t) + y1(t-5)
+%# d y3(t)/dt  = -y3(t) + y2(t-10)*y1(t-10)

+%# @end example
+%#
+%# one might do the following
+%#
+%# @example
+%# function f = fun (t, y, yd)

+%# f(1) =-y(1);                   %% y1' = -y1(t)
+%# f(2) =-y(2) + yd(1,1);         %% y2' = -y2(t) + y1(t-lags(1))
+%# f(3) =-y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) + y2(t-lags(2))*y1(t-lags(2))
+%# endfunction
+%# T = [0,20]

+%# res = ode54d (@fun, t, [1;1;1], [5, 10], ones (3,2));
+%# @end example
+%#
+%# @end itemize
 %# @end deftypefn
 %#
 %# @seealso{odepkg}
Index: ode78d.m
===================================================================
--- ode78d.m    (revision 6438)
+++ ode78d.m    (revision 6439)
@@ -23,11 +23,21 @@
 %#
 %# If this function is called with no return argument then plot the solution over time in a figure window while solving the set of DDEs that are defined in a function and specified by the function handle @var{@@fun}. The second input argument @var{slot} is a double vector that defines the time slot, @var{init} is a double vector that defines the initial values of the states, @var{lags} is a double vector that describes the lags of time, @var{hist} is a double matrix and describes the history of the DDEs, @var{opt} can optionally be a structure array that keeps the options created with the command @command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be other input arguments of any type that have to be passed to the function defined by @var{@@fun}.
 %#
+%# In other words, this function will solve a problem of the form
+%# @example
+%# dy/dt = fun (t, y(t), y(t-lags(1), y(t-lags(2), @dots{})))
+%# y(slot(1)) = init
+%# y(slot(1)-lags(1)) = hist(1), y(slot(1)-lags(2)) = hist(2), @dots{}
+%# @end example
+%#
 %# If this function is called with one return argument then return the solution @var{sol} of type structure array after solving the set of DDEs. The solution @var{sol} has the fields @var{x} of type double column vector for the steps chosen by the solver, @var{y} of type double column vector for the solutions at each time step of @var{x}, @var{solver} of type string for the solver name and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector that keep the informations of the event function if an event function handle is set in the option argument @var{opt}.
 %#
 %# If this function is called with more than one return argument then return the time stamps @var{t}, the solution values @var{y} and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector.
 %#
-%# For example, solve an anonymous implementation of a chaotic behavior
+%# For example:
+%# @itemize @minus
+%# @item
+%# the following code solves an anonymous implementation of a chaotic behavior
 %# @example
 %# fcao = @@(vt, vy, vz) [2 * vz / (1 + vz^9.65) - vy];
 %#
@@ -37,6 +47,29 @@
 %# vlag = interp1 (vsol.x, vsol.y, vsol.x - 2);
 %# plot (vsol.y, vlag); legend ("fcao (t,y,z)");
 %# @end example
+%#
+%# @item
+%# to solve the following problem with two delayed state variables
+%#
+%# @example
+%# d y1(t)/ dt = -y1(t)

+%# d y2(t)/ dt = -y2(t) + y1(t-5)
+%# d y3(t)/dt  = -y3(t) + y2(t-10)*y1(t-10)

+%# @end example
+%#
+%# one might do the following
+%#
+%# @example
+%# function f = fun (t, y, yd)

+%# f(1) =-y(1);                   %% y1' = -y1(t)
+%# f(2) =-y(2) + yd(1,1);         %% y2' = -y2(t) + y1(t-lags(1))
+%# f(3) =-y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) + y2(t-lags(2))*y1(t-lags(2))
+%# endfunction
+%# T = [0,20]

+%# res = ode78d (@fun, t, [1;1;1], [5, 10], ones (3,2));
+%# @end example
+%#
+%# @end itemize
 %# @end deftypefn
 %#
 %# @seealso{odepkg}
Index: ode23d.m
===================================================================
--- ode23d.m    (revision 6438)
+++ ode23d.m    (revision 6439)
@@ -23,11 +23,22 @@
 %#
 %# If this function is called with no return argument then plot the solution over time in a figure window while solving the set of DDEs that are defined in a function and specified by the function handle @var{@@fun}. The second input argument @var{slot} is a double vector that defines the time slot, @var{init} is a double vector that defines the initial values of the states, @var{lags} is a double vector that describes the lags of time, @var{hist} is a double matrix and describes the history of the DDEs, @var{opt} can optionally be a structure array that keeps the options created with the command @command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be other input arguments of any type that have to be passed to the function defined by @var{@@fun}.
 %#
+%# In other words, this function will solve a problem of the form
+%# @example
+%# dy/dt = fun (t, y(t), y(t-lags(1), y(t-lags(2), @dots{})))
+%# y(slot(1)) = init
+%# y(slot(1)-lags(1)) = hist(1), y(slot(1)-lags(2)) = hist(2), @dots{}
+%# @end example
+%#
 %# If this function is called with one return argument then return the solution @var{sol} of type structure array after solving the set of DDEs. The solution @var{sol} has the fields @var{x} of type double column vector for the steps chosen by the solver, @var{y} of type double column vector for the solutions at each time step of @var{x}, @var{solver} of type string for the solver name and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector that keep the informations of the event function if an event function handle is set in the option argument @var{opt}.
 %#
 %# If this function is called with more than one return argument then return the time stamps @var{t}, the solution values @var{y} and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector.
 %#
-%# For example, solve an anonymous implementation of a chaotic behavior
+%# For example:
+%# @itemize @minus
+%# @item
+%# the following code solves an anonymous implementation of a chaotic behavior
+%#
 %# @example
 %# fcao = @@(vt, vy, vz) [2 * vz / (1 + vz^9.65) - vy];
 %#
@@ -37,6 +48,29 @@
 %# vlag = interp1 (vsol.x, vsol.y, vsol.x - 2);
 %# plot (vsol.y, vlag); legend ("fcao (t,y,z)");
 %# @end example
+%#
+%# @item
+%# to solve the following problem with two delayed state variables
+%#
+%# @example
+%# d y1(t)/ dt = -y1(t)

+%# d y2(t)/ dt = -y2(t) + y1(t-5)
+%# d y3(t)/dt  = -y3(t) + y2(t-10)*y1(t-10)

+%# @end example
+%#
+%# one might do the following
+%#
+%# @example
+%# function f = fun (t, y, yd)

+%# f(1) =-y(1);                   %% y1' = -y1(t)
+%# f(2) =-y(2) + yd(1,1);         %% y2' = -y2(t) + y1(t-lags(1))
+%# f(3) =-y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) + y2(t-lags(2))*y1(t-lags(2))
+%# endfunction
+%# T = [0,20]

+%# res = ode23d (@fun, t, [1;1;1], [5, 10], ones (3,2));
+%# @end example
+%#
+%# @end itemize
 %# @end deftypefn
 %#
 %# @seealso{odepkg}

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

On 7 Nov 2009, at 19:13, franco basaglia wrote:

> Thank you all for your work.
>
> I have a last question.
> I'd convert my ode system to a discrete-time system
> with annual time step. In thsi way I can simplify my delay question.
>
>  What is the octave solver to use? How can I plan the system in  
> Octave?
>
> The discrete system is something like that:
>
> y1(t) = y1(t-1) -y1(t)
> y2(t) = y2(t-1) -y2(t) + y1(t-5)
> y3(t) = t3(t-1) -y3(t) + y2(t-10)*y1(t-10)
>

what you have written above is the Backward-Euler
method applied to your differential system.
to implement it in Octave, you need just to:

1) initialize your state variables for t<=10
2) add a cycle over t=11:T
3) solve the non-linear system at each step with e.g. "fsolve"

function res = fun (y,ym1,ym5,ym10)
res(1) = -y(1) + ym1(1) -y(1);
res(2) = -y(2) + ym1(2) -y(2) + ym5(1);
res(3) = -y(3) + ym1(3) -y(3) + ym10(2)*ym10(1);
endfunction

y = ones(3,10);
for t = 11:100
y(:,t) = fsolve(@(x) fun(x,y(:,t-1),y(:,t-5),y(:,t-10)), y(:,t-1));
endfor

and you're all set.
anyway this is more or less equivalent to

function f = fun (t, y, yd)
f(1) =-y(1);
f(2) =-y(2) + yd(1,1);
f(3) =-y(3) + yd(2,2)*yd(1,2);
endfunction
T = [0,20];
res = ode45d (@fun, T, [1;1;1], [5, 10], ones (3,2));

but much less accurate, as you can see by running both implementations  
above and then typing:

plot(-10:89, y)
hold on
plot (res.x, res.y, 'x-')

it seems to me that, while it more or less captures the asymtotic  
behaviuor,
your simplified solver damps away the initial oscillations completely.

> Thank you
> best regards
> f.b.

HTH,
c.
_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[Carlo de Falco-2]

On 7 Nov 2009, at 20:08, Carlo de Falco wrote:

>
> On 7 Nov 2009, at 19:13, franco basaglia wrote:
>
>> Thank you all for your work.
>>
>> I have a last question.
>> I'd convert my ode system to a discrete-time system
>> with annual time step. In thsi way I can simplify my delay question.
>>
>> What is the octave solver to use? How can I plan the system in  
>> Octave?
>>
>> The discrete system is something like that:
>>
>> y1(t) = y1(t-1) -y1(t)
>> y2(t) = y2(t-1) -y2(t) + y1(t-5)
>> y3(t) = t3(t-1) -y3(t) + y2(t-10)*y1(t-10)
>>
>
> what you have written above is the Backward-Euler
> method applied to your differential system.
> to implement it in Octave, you need just to:
>
> 1) initialize your state variables for t<=10
> 2) add a cycle over t=11:T
> 3) solve the non-linear system at each step with e.g. "fsolve"
>
> function res = fun (y,ym1,ym5,ym10)
> res(1) = -y(1) + ym1(1) -y(1);
> res(2) = -y(2) + ym1(2) -y(2) + ym5(1);
> res(3) = -y(3) + ym1(3) -y(3) + ym10(2)*ym10(1);
> endfunction
>
> y = ones(3,10);
> for t = 11:100
> y(:,t) = fsolve(@(x) fun(x,y(:,t-1),y(:,t-5),y(:,t-10)), y(:,t-1));
> endfor
>
> and you're all set.
> anyway this is more or less equivalent to
>
> function f = fun (t, y, yd)
> f(1) =-y(1);
> f(2) =-y(2) + yd(1,1);
> f(3) =-y(3) + yd(2,2)*yd(1,2);
> endfunction
> T = [0,20];
> res = ode45d (@fun, T, [1;1;1], [5, 10], ones (3,2));
>
> but much less accurate, as you can see by running both  
> implementations above and then typing:
>
> plot(-10:89, y)
> hold on
> plot (res.x, res.y, 'x-')

sorry, I guess

 >> plot(-9:90, y)
 >> hold on
 >> plot (res.x, res.y, 'x-')

is more correct

> it seems to me that, while it more or less captures the asymtotic  
> behaviuor,
> your simplified solver damps away the initial oscillations completely.
>
>> Thank you
>> best regards
>> f.b.
>
> HTH,
> c.

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave

[franco basaglia]

Thanks a lot Carlo.

The is the best way.
Using  discrete system I have exactly the results like Stella.
I'm happy


best regards
f.b.

2009/11/7 Carlo de Falco <carlo.defalco@...>

On 7 Nov 2009, at 19:13, franco basaglia wrote:

Thank you all for your work.

I have a last question.
I'd convert my ode system to a discrete-time system
with annual time step. In thsi way I can simplify my delay question.

 What is the octave solver to use? How can I plan the system in Octave?

The discrete system is something like that:

y1(t) = y1(t-1) -y1(t)
y2(t) = y2(t-1) -y2(t) + y1(t-5)
y3(t) = t3(t-1) -y3(t) + y2(t-10)*y1(t-10)

what you have written above is the Backward-Euler
method applied to your differential system.
to implement it in Octave, you need just to:

1) initialize your state variables for t<=10
2) add a cycle over t=11:T
3) solve the non-linear system at each step with e.g. "fsolve"

function res = fun (y,ym1,ym5,ym10)
res(1) = -y(1) + ym1(1) -y(1);
res(2) = -y(2) + ym1(2) -y(2) + ym5(1);
res(3) = -y(3) + ym1(3) -y(3) + ym10(2)*ym10(1);
endfunction

y = ones(3,10);
for t = 11:100
y(:,t) = fsolve(@(x) fun(x,y(:,t-1),y(:,t-5),y(:,t-10)), y(:,t-1));
endfor

and you're all set.
anyway this is more or less equivalent to

function f = fun (t, y, yd)
f(1) =-y(1);

f(2) =-y(2) + yd(1,1);

f(3) =-y(3) + yd(2,2)*yd(1,2);

endfunction
T = [0,20];
res = ode45d (@fun, T, [1;1;1], [5, 10], ones (3,2));

but much less accurate, as you can see by running both implementations above and then typing:

plot(-10:89, y)
hold on
plot (res.x, res.y, 'x-')

it seems to me that, while it more or less captures the asymtotic behaviuor,
your simplified solver damps away the initial oscillations completely.

Thank you
best regards
f.b.

HTH,
c.

_______________________________________________
Help-octave mailing list
Help-octave@...
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave