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