You can remove the error message by resetting the integrator at the start of the loop. Just do this:

tf := 100:
dsn(1e-16):
interface(warnlevel=0):
A := Array([0]):
while A[-1] < tf do
  dsn(tf):
  A(ArrayNumElems(A)+1) := dsn(eventfired=[1])[1]:
  dsn(eventclear):
  end do:
interface(warnlevel=3):    
A := convert(A,list):

The call 'dsn(1e-16)' resets the integrator so that it knows it will be integrating starting at zero to the right, so this also takes care of 'direction' as well.

- Allan

You can remove the error message by resetting the integrator at the start of the loop. Just do this:

tf := 100:
dsn(1e-16):
interface(warnlevel=0):
A := Array([0]):
while A[-1] < tf do
  dsn(tf):
  A(ArrayNumElems(A)+1) := dsn(eventfired=[1])[1]:
  dsn(eventclear):
  end do:
interface(warnlevel=3):    
A := convert(A,list):

The call 'dsn(1e-16)' resets the integrator so that it knows it will be integrating starting at zero to the right, so this also takes care of 'direction' as well.

- Allan

Hello Patrick,

If you need to get every time at which the event has fired, then you need to use halt events in combination with event cleearing.

The bolded comment refers to returning the time at which all events last fired - the list is for multiple events. The internal state machine only stores the last eveent fired point. We have many problems where events fire millions of times, so as a general rule storing all the fired points can pretty much lock up your computer from the storage requirement.

So first, here's how direction works:

> dsn := dsolve({diff(y(t),t)=-1, y(0)=1, s(0)=0}, numeric,
>    discrete_variables=[s(t)], events=[
>       [y(t), [y(t)=1, halt]]
>    ]);
                      dsn := proc(x_rkf45)  ...  end proc

> dsn(eventfired=[1]);
Error, (in _dat[1]) 'direction' must be set prior to calling/setting
'eventfired'
> dsn(direction=1):
> dsn(eventfired=[1]);
                                     [0.]

Now say you want to acquire the sequence of times at which the events fire. Here's how:

> dsn(1e10);
Warning, cannot evaluate the solution further right of 1.0000000, event #1
triggered a halt
                 [t = 1.00000000000000, y(t) = 1., s(t) = 0.]

> pts := dsn(eventfired=[1])[1];
                            pts := 1.00000000000000

> dsn(eventclear);
> dsn(1e10);                   
Warning, cannot evaluate the solution further right of 2.0000000, event #1
triggered a halt
                 [t = 2.00000000000000, y(t) = 1., s(t) = 0.]

> pts := pts,dsn(eventfired=[1])[1];
                   pts := 1.00000000000000, 2.00000000000000

> dsn(eventclear);                 
> dsn(1e10);                       
Warning, cannot evaluate the solution further right of 3.0000000, event #1
triggered a halt
                 [t = 3.00000000000000, y(t) = 1., s(t) = 0.]

> pts := pts,dsn(eventfired=[1])[1];
          pts := 1.00000000000000, 2.00000000000000, 3.00000000000000

> dsn(eventclear);                 
> dsn(1e10);                       
Warning, cannot evaluate the solution further right of 4.0000000, event #1
triggered a halt
                 [t = 4.00000000000000, y(t) = 1., s(t) = 0.]

> pts := pts,dsn(eventfired=[1])[1];
 pts := 1.00000000000000, 2.00000000000000, 3.00000000000000, 4.00000000000000

So you set the event to a halt event, store the value at the halt, then clear the event, rerun, repeat.

- Allan

Hello Patrick,

If you need to get every time at which the event has fired, then you need to use halt events in combination with event cleearing.

The bolded comment refers to returning the time at which all events last fired - the list is for multiple events. The internal state machine only stores the last eveent fired point. We have many problems where events fire millions of times, so as a general rule storing all the fired points can pretty much lock up your computer from the storage requirement.

So first, here's how direction works:

> dsn := dsolve({diff(y(t),t)=-1, y(0)=1, s(0)=0}, numeric,
>    discrete_variables=[s(t)], events=[
>       [y(t), [y(t)=1, halt]]
>    ]);
                      dsn := proc(x_rkf45)  ...  end proc

> dsn(eventfired=[1]);
Error, (in _dat[1]) 'direction' must be set prior to calling/setting
'eventfired'
> dsn(direction=1):
> dsn(eventfired=[1]);
                                     [0.]

Now say you want to acquire the sequence of times at which the events fire. Here's how:

> dsn(1e10);
Warning, cannot evaluate the solution further right of 1.0000000, event #1
triggered a halt
                 [t = 1.00000000000000, y(t) = 1., s(t) = 0.]

> pts := dsn(eventfired=[1])[1];
                            pts := 1.00000000000000

> dsn(eventclear);
> dsn(1e10);                   
Warning, cannot evaluate the solution further right of 2.0000000, event #1
triggered a halt
                 [t = 2.00000000000000, y(t) = 1., s(t) = 0.]

> pts := pts,dsn(eventfired=[1])[1];
                   pts := 1.00000000000000, 2.00000000000000

> dsn(eventclear);                 
> dsn(1e10);                       
Warning, cannot evaluate the solution further right of 3.0000000, event #1
triggered a halt
                 [t = 3.00000000000000, y(t) = 1., s(t) = 0.]

> pts := pts,dsn(eventfired=[1])[1];
          pts := 1.00000000000000, 2.00000000000000, 3.00000000000000

> dsn(eventclear);                 
> dsn(1e10);                       
Warning, cannot evaluate the solution further right of 4.0000000, event #1
triggered a halt
                 [t = 4.00000000000000, y(t) = 1., s(t) = 0.]

> pts := pts,dsn(eventfired=[1])[1];
 pts := 1.00000000000000, 2.00000000000000, 3.00000000000000, 4.00000000000000

So you set the event to a halt event, store the value at the halt, then clear the event, rerun, repeat.

- Allan

Hello Patrick,

The tobegin and toend directives are specialized to event iteration, which is specific to discrete equation state machine event iteration. This is a very specialized (niche) aspect of the tool, and based on the type of usage I have seen, you don't need it.

You want to use 'eventfired', which is described on the page, but lacking an example:

> dsn := dsolve({diff(y(t),t)=-1, y(0)=1, s(0)=0}, numeric,
>    discrete_variables=[s(t)], events=[
>       [y(t), [y(t)=1, s(t)=s(t)+y(t)]]
>    ]);
                      dsn := proc(x_rkf45)  ...  end proc

> dsn(4.5);
                [t = 4.5, y(t) = 0.500000000000001, s(t) = 4.]

> dsn(eventfired=[1]);
                              [4.00000000000000]

> dsn(6.1);
                [t = 6.1, y(t) = 0.900000000000001, s(t) = 6.]

> dsn(eventfired=[1]);
                              [6.00000000000000]

Using the query form of eventfired allows one to ask when any event last fired.

To query multiple events, say events 1,5,12, simply call:

> dsn(eventfired=[1,5,12]);

and the times will be provided.

If the event has not fired, you get the initial point.

- Allan

Hello Patrick,

The tobegin and toend directives are specialized to event iteration, which is specific to discrete equation state machine event iteration. This is a very specialized (niche) aspect of the tool, and based on the type of usage I have seen, you don't need it.

You want to use 'eventfired', which is described on the page, but lacking an example:

> dsn := dsolve({diff(y(t),t)=-1, y(0)=1, s(0)=0}, numeric,
>    discrete_variables=[s(t)], events=[
>       [y(t), [y(t)=1, s(t)=s(t)+y(t)]]
>    ]);
                      dsn := proc(x_rkf45)  ...  end proc

> dsn(4.5);
                [t = 4.5, y(t) = 0.500000000000001, s(t) = 4.]

> dsn(eventfired=[1]);
                              [4.00000000000000]

> dsn(6.1);
                [t = 6.1, y(t) = 0.900000000000001, s(t) = 6.]

> dsn(eventfired=[1]);
                              [6.00000000000000]

Using the query form of eventfired allows one to ask when any event last fired.

To query multiple events, say events 1,5,12, simply call:

> dsn(eventfired=[1,5,12]);

and the times will be provided.

If the event has not fired, you get the initial point.

- Allan

Hi Patrick,

The part that follows the trigger is a list of assignments that are to be performed when the trigger fires.

So when the trigger first fires at t=1 and y(1)=0, we get:

y(1) := 1:

s1(1) := s1(1)+y(1) = 0+1 = 0

s2(1) := s2(1)+pre(y(1)) = 0+0 = 0

Now since y is reset to 1, another event will fire at t=2, which will do:

y(2) := 1:

s1(2) := s1(2)+y(2) = 1+1 = 2

s2(2) := s2(2)+pre(y(2)) = 0+0 = 0

and so on.

If one were to plot the value of y(t), it would look like a sawtooth. The value of s1(t) is identically zero, and the value of s2(t) looks like floor(t).

These sequences of assignments allow one to program certain actions to occur during events.

- Allan

Hi Patrick,

The part that follows the trigger is a list of assignments that are to be performed when the trigger fires.

So when the trigger first fires at t=1 and y(1)=0, we get:

y(1) := 1:

s1(1) := s1(1)+y(1) = 0+1 = 0

s2(1) := s2(1)+pre(y(1)) = 0+0 = 0

Now since y is reset to 1, another event will fire at t=2, which will do:

y(2) := 1:

s1(2) := s1(2)+y(2) = 1+1 = 2

s2(2) := s2(2)+pre(y(2)) = 0+0 = 0

and so on.

If one were to plot the value of y(t), it would look like a sawtooth. The value of s1(t) is identically zero, and the value of s2(t) looks like floor(t).

These sequences of assignments allow one to program certain actions to occur during events.

- Allan

Hi Patrick,

It is possible to define a sort-of state machine that runs when an event occurs. Essentially a group of discrete equations that activate when any event fires. This is where the use of pre is relevant. The state machine can change the value of 'y(t)' on the fly, so the 'pre' allows an unambiguous reference to a changing variable.

Here is a simple example:

> dsn := dsolve({diff(y(t),t)=-1, y(0)=1, s1(0)=0, s2(0)=0}, numeric,
>    discrete_variables=[s1(t),s2(t)], events=[
>       [y(t), [y(t)=1, s1(t)=s1(t)+y(t), s2(t)=s2(t)+pre(y(t))]]
>    ]);
                      dsn := proc(x_rkf45)  ...  end proc

> dsn(5.5);                                                         
[t = 5.5, y(t) = 0.500000000000001, s1(t) = 5.,

                                 -15
    s2(t) = -0.938946186426759 10   ]

So this shows that the sum s1 without the pre is getting the updated value of y, while the sum s2 with the pre is getting the value of y at the beginning of the event.

- Allan

Hi Patrick,

It is possible to define a sort-of state machine that runs when an event occurs. Essentially a group of discrete equations that activate when any event fires. This is where the use of pre is relevant. The state machine can change the value of 'y(t)' on the fly, so the 'pre' allows an unambiguous reference to a changing variable.

Here is a simple example:

> dsn := dsolve({diff(y(t),t)=-1, y(0)=1, s1(0)=0, s2(0)=0}, numeric,
>    discrete_variables=[s1(t),s2(t)], events=[
>       [y(t), [y(t)=1, s1(t)=s1(t)+y(t), s2(t)=s2(t)+pre(y(t))]]
>    ]);
                      dsn := proc(x_rkf45)  ...  end proc

> dsn(5.5);                                                         
[t = 5.5, y(t) = 0.500000000000001, s1(t) = 5.,

                                 -15
    s2(t) = -0.938946186426759 10   ]

So this shows that the sum s1 without the pre is getting the updated value of y, while the sum s2 with the pre is getting the value of y at the beginning of the event.

- Allan

The problem with rkf45 is that it is a lower order method, and a tolerance of 1e-12 is a little tight for the working precision, so the solution may fail for some probems. Loosening it to 1e-11 should work fine.

There is one additional detail that I did not think about when sending you the modified procedure, as I did not consider the need for use of dverk78 for the first solution. You can provide a much more efficient solution by combining the original problem (sys_Ode) with the other problems one at a time, forming several larger systems that can be solved simultaneously with the same method.

This approach will have greater accuracy, as any large effects that a smooth change in the original solution will have on the derived solution (the P0, Omega, ...) will automatically result in a step reduction for both solutions (as they are being simult. solved).

Here is the code using only dverk78 for all solutions:

Digits:=15;

gradPhi4:= proc (P0,Gamma,Omega,lambda)
local A,B,ics,sys,sol,Pv,
   sysP0,PP0,LP0,WP0,dsolP0,P_P0,
   sysGamma,PGamma,LGamma,WGamma,dsolGamma,P_Gamma,
   sysOmega,POmega,LOmega,WOmega,dsolOmega,P_Omega,
   syslambda,Plambda,Llambda,Wlambda,dsollambda,P_lambda,
   opts,Pop,time,a,b,c,d,t;

   opts := numeric, method=dverk78, relerr=1e-10, abserr=1e-10,
      output=listprocedure;

   A := 1e-7; B:= 1e4;
   sys := {
      diff(P(t),t)= A*Gamma*P(t)*(L(t)-P(t))-A*Omega*W(t)*P(t),
      diff(L(t),t)= -P(t),
      diff(W(t),t)=P(t),
      P(0)=P0, L(0)=B*lambda, W(0)=0
   };
   sol := dsolve(sys, opts);
   Pv := rhs(sol[3]);

   sysP0 := {
      diff(PP0(t),t)=A*Gamma*(PP0(t)*(L(t)-P(t))+P(t)*(LP0(t)-PP0(t)))
                     -A*Omega*(WP0(t)*P(t)+W(t)*PP0(t)),
      diff(LP0(t),t)=-PP0(t),
      diff(WP0(t),t)=PP0(t),
      PP0(0)=1, LP0(0)=0, WP0(0)=0
   }:
   dsolP0 := dsolve(sysP0 union sys, opts);
   P_P0 := eval(PP0(t),dsolP0[2..-1]);

   sysGamma := {
      diff(PGamma(t),t)=A*P(t)*(L(t)-P(t))
         +A*Gamma*(PGamma(t)*(L(t)-P(t))+P(t)*(LGamma(t)-PGamma(t)))
         -A*Omega*(WGamma(t)*P(t)+W(t)*PGamma(t)),
      diff(LGamma(t),t)=-PGamma(t),
      diff(WGamma(t),t)=PGamma(t),
      PGamma(0)=0, LGamma(0)=0, WGamma(0)=0
   }:
   dsolGamma:= dsolve(sysGamma union sys, opts);
   P_Gamma := eval(PGamma(t),dsolGamma[2..-1]);

   sysOmega := {
      diff(POmega(t),t)=A*Gamma*(POmega(t)*(L(t)-P(t))+P(t)*(LOmega(t)
         -POmega(t)))-A*W(t)*P(t)-A*Omega*(WOmega(t)*P(t)+W(t)*POmega(t)),
      diff(LOmega(t),t)=-POmega(t),
      diff(WOmega(t),t)=POmega(t),
      POmega(0)=0, LOmega(0)=0, WOmega(0)=0
   }:
   dsolOmega := dsolve(sysOmega union sys, opts);
   P_Omega := eval(POmega(t),dsolOmega[2..-1]);

   syslambda := {
      diff(Plambda(t),t)=A*Gamma*(Plambda(t)*(L(t)-P(t))+P(t)*(Llambda(t)
         -Plambda(t)))-A*Omega*(Wlambda(t)*P(t)+W(t)*Plambda(t)),
      diff(Llambda(t),t)=-Plambda(t),
      diff(Wlambda(t),t)=Plambda(t),
      Plambda(0)=0, Llambda(0)=B, Wlambda(0)=0
   }:
   dsollambda := dsolve(syslambda union sys, opts);
   P_lambda:=eval(Plambda(t),dsollambda[2..-1]);

   Pop:=[0.333333,0.333333,0.166667,0.666667,0.166667,0.166667,1.833333,
      1.166667,1.833333,2.5,12.5,27.66667,60.66667,45,81.5,99.833333,23.66667,
      81.16667,59,48.83333];
   time:=[30,90,150,180,210,240,270,330,360,390,420,450,480,510,540,570,600,
      630,660,690];

   a:=add((Pv(time[i])-Pop[i])*P_P0(time[i]),i=1..20);
   b:=add((Pv(time[i])-Pop[i])*P_Gamma(time[i]),i=1..20);
   c:=add((Pv(time[i])-Pop[i])*P_Omega(time[i]),i=1..20);
   d:=add((Pv(time[i])-Pop[i])*P_lambda(time[i]),i=1..20);

   Vector[column](1..4,[a,b,c,d]);
end proc:

And for this code, I get the result:

   [-0.474777084408743]
   [-0.996172632304337]
   [-0.500013706760093]
   [-0.990140063780814]

Note that I had to loosen the tolerances to 1e-10 for this to succeed, which hints to me that there is some non-trivial interaction between the two subsystems that is likely being missed by the 2-step approach.

Also, if you know in advance that the solution is smooth (which seems likely as the inputs are all polynomial), then you may be able to get even higher accuracy relatively efficiently using Maple's gear method - this or taylorseries are the methods I use when I want to compute high accuacy solutions for ODE problems.

The gear method is a variable order Richardson-extrapolation related technique, and can give very accurate results.

Using gear I get:

15 Digits: abs/relerr = 1e-13 (0.7 sec)
   [-0.47477705310  ]
   [-0.996172630850 ]
   [-0.5000137068574]        
   [-0.990140058009 ]

20 Digits: abs/relerr = 1e-16 (20 sec)
   [-0.4747770821795167  ]
   [-0.99617263220360641 ]
   [-0.500013706766194511]
   [-0.99014006337411834 ]

25 Digits: abs/relerr = 1e-20 (67 sec)
   [-0.474777082169271469957  ]
   [-0.9961726322030855172671 ]
   [-0.50001370676622364282772]
   [-0.9901400633724119371884 ]           
 

which suggests that the 20 Digit result is accurate to about 11 Digits, and obtained in only 20 sec., so one might expect roughly 15 Digits accuracy in the last result.

The dverk78 method is a little different than the rkf45 method, as it uses a continuous error control scheme that bounds the error throughout the integration step using a continous interpolant. This type of error control requires significantly more function evaluations than a straight up rk78 method, but it does an excellent job of detecting prroblem points in the domain of integration. The addition of the piecewise and range options was only done to for the two main default solvers in Maple, the non-stiff default (rkf45) and the stiff default (rosenbrock). The addition of these options to dverk78 would be no small task, as the solver was not written with solution storage in mind (while rkf45 and rosenbrock were).

As for why you need to set both abserr and relerr, this is explained in the help page ?dsolve,Error_Control, but the gist of it is that we accept an integration step when:

error_estimate < abserr + relerr*|solution|

so only reducing the abserr allows for a relative error of size relerr, and similarly only reducing the relerr allows for an absolute error of abserr. In order to really reduce the error in all cases, both control values need to be scaled down.

Hope this answers your questions!

Allan Wittkopf

Math Developer, Maplesoft

Now that I have had a better chance to look at the procedure, and the problem, I can better explain what you are seeing here.

First of all, the values you are comparing between Maple 11 and Maple 13 are not the values of the ODE solution, but rather differences between the ODE solution values and a nearby exact/measured solution.

If the exact/measured solution is close to the values that should be produced by the numerical solution, then small changes in the numerical solution will give large changes to the difference.

As a simple illustration, the number 3.14159 is Pi computed accurately to 6 Digits (read as a relative error tolerance of 1e-6), but then compare 3.14159-Pi. The latter difference has magnitude ~ 2.65e-6, but a 1e-6 relative difference in the computed value from 3.14159 to 3.14160 makes an O(1) change to the difference 3.14160-Pi to 7.35e-6.

The problem you showed is running the 4 main integrations with a relative error tolerance bound of 1e-6, but the value being computed is a difference between the solution and nearby values, so this amplifies the difference.

Now as to why there is a difference - the error control mechanism for rkf45 was made more stringent for Maple 13 than it was in prior versions - there were some classes of problems where the old error control was insufficient. This means that in Maple 13, the solutions are being computed with greater accuracy than they were in Maple 11.

To demonstrate that this is true, I re-ran your example code adding in the option 'relerr=1e-7' to the 4 rkf solver calls in the code, and the overall results are below:

Maple 11 original:

   [ -0.500658039487016757813 ]
   [-1.0000141405193850819475 ]
   [-0.49999905745147104872066]
   [-1.0009082015824584458607 ]

Maple 11 w/relerr=1e-7

   [ -0.475420032178772600936 ]
   [-0.9961429623618422930957 ]
   [-0.50000737987411704927995]
   [-0.9915922526684901284710 ]

Maple13:

   [ -0.476619191679594384338 ]
   [-0.9965175018113306636274 ]
   [-0.50008152041009398631850]
   [-0.9889420779552149907234 ]
 

The main things to note are that the Maple 11 solution changes significantly when the relative error is changed from 1e-6 (the default) to 1e-7, and that the Maple 11 result with 1e-7 is closer to the Maple 13 result.

Now in the prior email I made a few comments, so I modified your code to follow some of the ideas in the comments, namely:

- Compute the dverk78 solution using rkf45 and piecewise output instead, so that no recursion is needed for the subsequent 4 problems (efficiemcy)

- Tighten the relative error tolerances, and loosen the absolute (accuracy).

- Perform the computation in 15 Digits, allowing fast hardware computation to be used (efficiency).

With these changes the run-time of the problem drops from ~400 sec. to ~16 sec., and the accuracy of the solution appears to be roughly 5-6 Digits.

Here is the modified code:

Digits:=15;

gradPhi4:= proc (P0,Gamma,Omega,lambda)
local A,B,ics,sys_Ode,sol,P,L,W,
   sysP0,PP0,LP0,WP0,dsolP0,P_P0,
   sysGamma,PGamma,LGamma,WGamma,dsolGamma,P_Gamma,
   sysOmega,POmega,LOmega,WOmega,dsolOmega,P_Omega,
   syslambda,Plambda,Llambda,Wlambda,dsollambda,P_lambda,
   opts,Pop,time,a,b,c,d,t;

   A := 1e-7; B:= 1e4;
   ics := P(0)=P0, L(0)=B*lambda, W(0)=0;
   sys_Ode := diff(P(t),t)= A*Gamma*P(t)*(L(t)-P(t))-A*Omega*W(t)*P(t),
              diff(L(t),t)= -P(t),diff(W(t),t)=P(t);
   sol := dsolve([sys_Ode,ics], numeric, method=rkf45, abserr=1e-12,
      relerr=1e-12, output=piecewise, range=0..750);
   P := rhs(sol[3]); L := rhs(sol[2]); W := rhs(sol[4]);
   P := subsop(1=NULL,2=NULL,-1=NULL,P);
   L := subsop(1=NULL,2=NULL,-1=NULL,L);
   W := subsop(1=NULL,2=NULL,-1=NULL,W);

   opts := numeric, relerr=1e-12, abserr=1e-12, output=listprocedure;

   sysP0 := {
      diff(PP0(t),t)=A*Gamma*(PP0(t)*(L-P)+P*(LP0(t)-PP0(t)))
                     -A*Omega*(WP0(t)*P+W*PP0(t)),
      diff(LP0(t),t)=-PP0(t),
      diff(WP0(t),t)=PP0(t),
      PP0(0)=1, LP0(0)=0, WP0(0)=0
   }:
   dsolP0 := dsolve(sysP0, opts);
   P_P0 := eval(PP0(t),dsolP0[2..-1]);

   sysGamma := {
      diff(PGamma(t),t)=A*P*(L-P)
         +A*Gamma*(PGamma(t)*(L-P)+P*(LGamma(t)-PGamma(t)))
         -A*Omega*(WGamma(t)*P+W*PGamma(t)),
      diff(LGamma(t),t)=-PGamma(t),
      diff(WGamma(t),t)=PGamma(t),
      PGamma(0)=0, LGamma(0)=0, WGamma(0)=0
   }:
   dsolGamma:= dsolve(sysGamma, opts);
   P_Gamma := eval(PGamma(t),dsolGamma[2..-1]);

   sysOmega := {
      diff(POmega(t),t)=A*Gamma*(POmega(t)*(L-P)+P*(LOmega(t)-POmega(t)))
         -A*W*P-A*Omega*(WOmega(t)*P+W*POmega(t)),
      diff(LOmega(t),t)=-POmega(t),
      diff(WOmega(t),t)=POmega(t),
      POmega(0)=0, LOmega(0)=0, WOmega(0)=0
   }:
   dsolOmega := dsolve(sysOmega, opts);
   P_Omega := eval(POmega(t),dsolOmega[2..-1]);

   syslambda := {
      diff(Plambda(t),t)=A*Gamma*(Plambda(t)*(L-P)+P*(Llambda(t)-Plambda(t)))
         -A*Omega*(Wlambda(t)*P+W*Plambda(t)),
      diff(Llambda(t),t)=-Plambda(t),
      diff(Wlambda(t),t)=Plambda(t),
      Plambda(0)=0, Llambda(0)=B, Wlambda(0)=0
   }:
   dsollambda := dsolve(syslambda, opts);
   P_lambda:=eval(Plambda(t),dsollambda[2..-1]);

   Pop:=[0.333333,0.333333,0.166667,0.666667,0.166667,0.166667,1.833333,
      1.166667,1.833333,2.5,12.5,27.66667,60.66667,45,81.5,99.833333,23.66667,
      81.16667,59,48.83333];
   time:=[30,90,150,180,210,240,270,330,360,390,420,450,480,510,540,570,600,
      630,660,690];

   a:=add((eval(P,t=time[i])-Pop[i])*P_P0(time[i]),i=1..20);
   b:=add((eval(P,t=time[i])-Pop[i])*P_Gamma(time[i]),i=1..20);
   c:=add((eval(P,t=time[i])-Pop[i])*P_Omega(time[i]),i=1..20);
   d:=add((eval(P,t=time[i])-Pop[i])*P_lambda(time[i]),i=1..20);

   Vector[column](1..4,[a,b,c,d]);
end proc:

And here is the result in M13:

   [ -0.47477721234 ]
   [-0.996172636229 ]
   [-0.5000137071106]        
   [-0.990140076778 ]        

Sincerley.

Allan Wittkopf

Math Developer, Maplesoft

I'm currently running your example procedure via:

gradPhi4( 0.03628943868909384264412104, 7.662243093848743041881986,  8.035227714257298233205133, 2.024474547998493950062639);

but had a few comments/questions:

1) Why are you using dverk78 for the first solution, but then using rkf45 for the subsequent solutions?

2) I was curious about the use of 'abserr=1e-15', but no specification of relerr. This results in solving a problem where the absolute error is limited to 1e-15, but the relative error is set to 1e-6 (rkf45) and 1e-8 (dverk78), so it is almost like running a pure relative error test. Is this the intention?

3) The fact that one solution with errors on the order of 1e-8 is being used as an input for another solution process will have the effect of amplifying the error, so a higher precision should probably be used when one solution is being used to compute another. Has the problem been symbolically analyzed to assure that this is not the case?

4) The efficiency of the solution may be a bit impacted by the fact that an 'on-the-fly' solution is being used to compute another. Use of one of the solutions with a piecewise or stored output for the first problem might be beneficial.
 

I will continue to look at the problem.

Sincerely,

Allan Wittkopf

Math Developer, Maplesoft

XP and Vista probably use the same math library, and they definitely use the exact same Maple binaries, so the results have a higher probability of being the same than on, say, 2 differing versions of Linux (where the math libraries are generally a little bit different).

If you compared 32-bit Windows XP to 64-bit Windows, I'd expect you'll probably see more of a difference.

Sincerly,

Allan Wittkopf

Math Developer, Maplesoft

Do you know if your ODE system is stable? i.e. do small perterbations in the initial conditions result in large changes to the results?

If so, then you may even find that you can get different results on different versions of Linux.

As a simple illustration, if the ODE system in question contained a sin(t) term, the external math libraries on the platform do not give full precision accuracy to the full hardware precision, and this may make a very small difference in one step of the integration, which for an unstable system, integrated over a sufficiently long time period, can make the solution very different.

Though you are not using hardware precision, the same rules might apply, as a number might be rounded up on one platform/flavor, and rounded down on another, resulting in a small diffference that can be amplified by an unstable system.

Without seeing/analyzing the system it's difficult to know whether this behavior is expected (as a result of the nature of the system) or a problem in one of the methods.