Re: octave and delay ode's

[Carlo de Falco]

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.