#81 num_odesolve does not handle intervals with pi. Here is the fix.

[Marduk Bolaños]

The package numeric provides num_int and num_odesolve. The former can handle integrals like

num_int(sin(x),x=(0 .. pi));.

However the following fails with ***** (*interval* 0 pi) invalid as numeric interval:

depend(x,t);
num_odesolve(df(x,t)=x, x=0, t=(0 .. pi));

The difference is that numint evaluates the interval in rounded mode. In the following, I implemented this for rungekuttaeval:

symbolic procedure rungekuttaeval u;
     % interface function;
  begin scalar e,f,x,y,sx,sy,en,d,v,w,q,z;
    u := for each x in u collect reval x;
    u := accuracycontrol(u,20,6);

     % equations
    e :=car u; u :=cdr u;
    e := if eqcar(e,'list) then cdr e else {e};

     % dependent variable
    q :=car u; u :=cdr u;
    q := if eqcar(q,'list) then cdr q else {q};
    for each yy in q do
    <<if not eqcar(yy,'equal) or not idp cadr yy
     then typerr(yy,"expression `dep. variable=starting value'");
     sy:=caddr yy.sy; y := cadr yy.y;
    >>;
    sy:=reversip sy; y:=reversip y;
    if length sy neq length e then
      rederr "number of equations and variables don't correspond";

     % independent variable
    x :=car u; u :=cdr u;
+   if not eqcar(x,'equal) then typerr(x,"interval bounds");
+   if not !*rounded then setdmode('rounded,z:=!*rounded:=t);
+   w:=revalnuminterval(caddr x,'infinity);
+   if z then setdmode('rounded,z:=!*rounded:=nil);
    sx:=car w; en:=cadr w; x := cadr x;
     % convert expressions to explicit ODE system.
     q := for each yy in y collect {'df,yy,x};
     v:=cdr solveeval list('list. e, 'list . q);
     if cdr v then
       ruksyserr(nil,"cannot convert to explicit ODE system");
     f := for each d in q collect rukufindrhs(d,v);
    return rungekutta1(f,x,y,sx,en,sy);
end;

I declared the variable z and added/modified the lines beginning with +.

[Marduk Bolaños]

Further testing reveals very weird behavior. It only works for some integer multiples of pi, others fail with +++ Error attempt to take car of an atom: nil.

[Marduk Bolaños]

It turns out that rungekuttares was not written carefully. Here is the corrected version, that checks the size of the list before applying car and cdr.

symbolic procedure rungekuttares(l,st);
   % eliminate intermediate points.
     if st=iterations!* then l else
     << for i:=1:iterations!* collect
         <<for j:=1:m do if (length l > 1) then l:= cdr l; car l>>
     >> where m=st/iterations!*;

[Marduk Bolaños]

By the way, the default accuracy, defined in the beginning of rungekuttaeval

u := accuracycontrol(u,20,6);

is larger than machine precision and very likely unnecessarily large.

My experiments with the harmonic oscillator show that with low accuracy (e.g. 6) after 500 oscillations and using 14,000 steps the energy is conserved. Which might also suggest that this method is sympletic. If someone knows the algorithm behind the code it would be good to explicitly mention it in the documentation.

[Rainer Schöpf]

I committed a correction for the the problem with pi as part of the interval specification.

I'm afraid that I cannot reproduce the car of nil error you mention. I suspect this is caused by rounding. Can you please supply a complete example?

Rainer

[Marduk Bolaños]

Now I can't reproduce it either. It is very strange.

[Rainer Schöpf]

Correction committed.