Allan Wittkopf

240 Reputation

7 Badges

18 years, 312 days

MaplePrimes Activity

These are answers submitted by Allan Wittkopf

First, regarding use of O(h) BCs with an O(h^2) method, this is viable because the O(h^2) method is being applied at O(1/h) points, so the overall error is O(h), and the BCs are at a finite number of points, so O(h) BCs will contribute O(h) to the overall error. That is the origin of the statement in the help.


Next item of note, the default scheme in use is the 1/2 theta method, or the Crank Nicholson method. The theta version is a parametrized version that gives:

theta=0 - Forward Euler, O(h)

theta=1/2 - Crank-Nicholson, O(h^2)

theta=1 - Backward Euler, O(h)

The higher order method was chosen as the default for pdsolve/numeric, rather than the more stable theta=1 method.

If you re-run your problem using 'method=theta[1]' you will notice that the numeric artifacts are gone.


The heat equation is a parabolic equation that has certain restrictions on the time/space steps to produce a fully stable solution with the O(h^2) scheme. I don't recall off the top of my head the exact conditions, but they are undesirable - something like requiring that the time step be < h^2/2, where h is the space step. It is far more reasonable to simply use the theta[1] method to achieve greater stability.

First it should be noted that in Maple, when an ODE solution is obtained, the values that the solution will be requested at (by default) are not known in advance, so the only possible method of operation is to compute solution values *as* they are requested.

The lsode method in Maple is more-or-less a direct mapping of the published Fortran algorithm with a Maple interface.

That interface was designed to integrate from a starting point to a requested ending point, outputting the value at that endpoint. The issue is that continuation of the integration (i.e. requesting of a second point) is performed from:

1) If the new point is farther away from the initial value than the last requested point, integration is performed from the last requested point

2) If the new point is closer to the initial value than the last requested point, integration is performed from the initial point.

One can actually observe this directly in Maple:

> mysin := proc(x) if not type(x,'numeric') then 'procname'(x);
>    else lprint(`Request at x=`,x); sin(x); end if; end proc:
> dsn := dsolve({diff(y(x),x)=mysin(x),y(0)=0}, numeric,
>    method=lsode, known={mysin});
`Request at x=`, 0
                      dsn := proc(x_lsode)  ...  end proc

> dsn(1);
`Request at x=`, 0.
`Request at x=`, .316227766016837939e-3
`Request at x=`, .316227766016837939e-3
`Request at x=`, .822526225536629685
`Request at x=`, .919835463478538284
`Request at x=`, 1.01714470142044688
                      [x = 1., y(x) = 0.459697513613591]

> dsn(1.2);
`Request at x=`, 1.
`Request at x=`, 1.00031207948782574
`Request at x=`, 1.00031207948782574
`Request at x=`, 1.13710328654956050
`Request at x=`, 1.19166312131692576
`Request at x=`, 1.24622295608429101
                      [x = 1.2, y(x) = 0.637642166314161]

> dsn(0.1);
`Request at x=`, 0.
`Request at x=`, .316227766016837953e-4
`Request at x=`, .632455532033675906e-4
`Request at x=`, .104069022625259522
`Request at x=`, .988159476246405072e-1
`Request at x=`, .119564316892035533
                     [x = 0.1, y(x) = 0.00499559764096088]

> dsn(1.2);
`Request at x=`, .100000000000000006
`Request at x=`, .100376805608346789
`Request at x=`, .100376805608346789
`Request at x=`, 1.03776125115700646
`Request at x=`, 1.13534129741291556
`Request at x=`, 1.28009936449751138
                      [x = 1.2, y(x) = 0.637642154095732]

Item 1 has the unfortunate side effect that the value computed for a point then depends upon all other points that have been requested before it, and can be changed if that sequence of points is different.

In fact this is directy apparent from the above example:

                      [x = 1.2, y(x) = 0.637642166314161]
                      [x = 1.2, y(x) = 0.637642154095732]

where the results at x=1.2 are different depending on whether we requested a value at 1.0 v.s. 0.1 before obtaining them. For some problems the differences can be much more substantial than the above.

This issue can be alleviated by the (brutal) use of the 'startinit=true' option (not really recommended) which forces consistency by performing a full integration from the initial point for each requested point.

In contrast, the default nonstiff (rkf45) and stiff (rosenbrock) methods are more recent Maple implementations with a tighter associtation with the Maple working style.

These methods are specifically implemented to resolve these issues (repeatability), as this was part of their design.

By default they use a natural interpolant associated with the integration method for point queries, and always take the 'natural' step (i.e. computed from the error criteria) to obtain a point request, so continuation past a previously computed solution value should always give the same result as computation dirrectly from the initial condition *WITHOUT* having to re-integrate from the initial point.

As for efficiency, there is a compile option that is supported by rkf45, rosenbrock and ck45 that is not available for the other solvers in Maple. The codes themselves are implemented in a compiled language, so the combination of the compiled integrators and the compilation of the problem should yield a hightly efficient solution process.

Here is a similar example as the above done with rkf45, where I only display the results of the calls:

                      [x = 1.0, y(x) = 0.459697678370392]
                      [x = 1.2, y(x) = 0.637642235621288]
                     [x = 0.1, y(x) = 0.00499585264319642]
                      [x = 1.2, y(x) = 0.637642235621288]

Note that the pair of results for 1.2 are *identical*.

Furthermore, the rkf45, rosenbrock, and ck45 methods *also* support a pre-compute/storage operating mode, where at the time of requesting the solution the integration is performed and the results are stored, making any subsequent solution value queries into a table lookup. This is controlled through use of the range option, but beware that sufficient storage must be allocated to accurately compute all solution values in the requested range using the requested tolerance, so memory usage may be a concern for a rapidly varying solution (or for use of a non-stiff method on a stiff problem).

To specify the recurrence in Maple you can enter it as follows:

> reqs := {
>    Y(n+1)=Z(n),
>    Z(n+1)=-alpha*(n+1)*Y(n)+(gamma*(n+1)+alpha*(n+1))*Z(n)
> };

Now Maple's recurrence solver is much stronger for solution of single recurrence equations rather than systems, especially for recurrences with non-constant coefficients, and in this case we can easily convert to this form as folows:

> req2 := subs(reqs[1],eval(reqs[2],n=n+1));
req2 :=

    Z(n + 2) = -alpha (n + 2) Z(n) + (gamma (n + 2) + alpha (n + 2)) Z(n + 1)

Where now we will specify Z(0), and Z(1)=Y(0) as the initial values for the recurrence.

For fixed values of the initial values and alpha and gamma, we can easily get a recurrence computation procedure as follows:

> rprc := rsolve(subs(alpha=-1/3,gamma=1\                                      
> /2,{req2}) union {Z(1)=1, Z(0)=1/2},Z(n),makeproc);
                        rprc := proc(n)  ...  end proc

> seq(rprc(i),i=0..10);
                                  100  196         3416  3094  80570
          1/2, 1, 2/3, 4/3, 16/9, ---, ---, 154/9, ----, ----, -----
                                  27   27           81    27    243

where we have obtained the first 11 terms of the recurrence, includng the 2 initial values.

For these values we can also obtain a closed form solution as follows:

> res := rsolve(subs(alpha=-1/3,gamma=1/\                                      
> 2,{req2}) union {Z(1)=1, Z(0)=1/2},Z(n));
            n                11        10          9            8
res := (1/6)  GAMMA(n + 2) (n   + 209 n   + 19140 n  + 1011450 n

                 7              6                5                 4
     + 34178133 n  + 773034801 n  + 11898687086 n  + 124103935780 n

                     3                  2                                    |
     + 855146138616 n  + 3683974520160 n  + 8859037402368 n + 8916100448256) |

                          35647153              10        9          8
    1/17832200896512 + -------------- (n + 2) (n   + 195 n  + 16422 n

               7             6              5               4                3
     + 783702 n  + 23371257 n  + 452814579 n  + 5738290928 n  + 46634083572 n

     + 230640572016 n  + 620505857664 n + 681091006464)

    /  n - 1              \
    | --------'           |
    |'  |  |    /    12  \|   /   11        10          9            8
    |   |  |    |- ------||  /  (n   + 209 n   + 19140 n  + 1011450 n
    |   |  |    \  _i + 3/| /
    |   |  |              |
    \  _i = 0             /

                 7              6                5                 4
     + 34178133 n  + 773034801 n  + 11898687086 n  + 124103935780 n

                     3                  2
     + 855146138616 n  + 3683974520160 n  + 8859037402368 n + 8916100448256)

                       /    -1               \
                       | --------'           |
          35647153     |'  |  |    /    12  \|
     - --------------- |   |  |    |- ------||
       445805022412800 |   |  |    \  _i + 3/|
                       |   |  |              |
                       \  _i = 0             /

       /n - 1  /               /  n1 - 1             \\\\
       |-----  |               | --------'           ||||
       | \     |   35647153    |'  |  |    /    12  \||||
     + |  )    |-------------- |   |  |    |- ------|||||
       | /     |68109100646400 |   |  |    \  _i + 3/||||
       |-----  |               |   |  |              ||||
       \n1 = 0 \               \  _i = 0             ////

where we can also evaluate the sums:

> res := value(res);
            n                11        10          9            8
res := (1/6)  GAMMA(n + 2) (n   + 209 n   + 19140 n  + 1011450 n

                 7              6                5                 4
     + 34178133 n  + 773034801 n  + 11898687086 n  + 124103935780 n

                     3                  2                                    |
     + 855146138616 n  + 3683974520160 n  + 8859037402368 n + 8916100448256) |

         18337          35647153              10        9          8
    - ------------ + -------------- (n + 2) (n   + 195 n  + 16422 n
      229323571200   34054550323200

               7             6              5               4                3
     + 783702 n  + 23371257 n  + 452814579 n  + 5738290928 n  + 46634083572 n

                     2                                       n   /    11
     + 230640572016 n  + 620505857664 n + 681091006464) (-12)   /  ((n

            10          9            8             7              6
     + 209 n   + 19140 n  + 1011450 n  + 34178133 n  + 773034801 n

                    5                 4                 3                  2
     + 11898687086 n  + 124103935780 n  + 855146138616 n  + 3683974520160 n

     + 8859037402368 n + 8916100448256) GAMMA(n + 3)) + 1/4903855246540800 |

                                   (n + 2)
    35647153 (1 + 11 exp(12)) (-12)

                       n                                           \
       5133190032 (-12)  (GAMMA(n + 3) - (n + 2) GAMMA(n + 2, -12))|   /
     - ------------------------------------------------------------|  /  (
                               GAMMA(n + 3)                        / /

         (n + 2)         |
    (-12)        exp(12))|

And we can then compare with the procedure form to verify the solution:

> eval(res,n=10),rprc(10);
                                 80570  80570
                                 -----, -----
                                  243    243

I tried to obtain a closed form solution for arbitrary values of alpha and gamma, but had little luck, but you may be able to obtain solutions for the values of those parameters that are of interest to you.

Hope this helps,


Even if the integrands in the call to dsolve did depend on 't' it is still possible to get a significant speedup for this problem.

There is a gotcha here that is avoided by the prior two answers, but it should be mentioned. When numeric dsolve calls the evaluation procedure, it does so with the hardware setting of Digits:

> evalhf(Digits);

The result of this is that the calls to evalf(Int()) are being asked to produce an answer that is accurate to 15 Digits, which means these calls need to compute integrand values with a setting of 17 or 18 Digits. This higher Digits setting puts Maple into software precision mode for the integrations, which is much slower than the hardware integrators.

The same timing can be achieved by reducing the requested precision to 10 Digits inside the procedure, which is what is used when the integrand values are computed outside of numeric dsolve anyway.

Note that the default tolerances for rkf45 are 1e-7, so a 10 Digit accurate integration is sufficient here.

In addition, use of multiple integration instead of recursive integration is recommended.

So the procedure simply needs to change to

dproc3 := proc (nn, t, Y, YP)
local n, nd, md;
for n to nn do
YP[n] := -((I*e^2*4)/(h*a^2))*Y[n]*add(add(Y[nd]*conjugate(Y[nd])
[y=0..a,x=0..a])),, nd=1..nn);
end do;
end proc:

Note though that if the integrands do not depend on 't', the prior replies provide the best approach.

Hello Patrick,

Sorry that the page is not easy to read - the problem is that there's just so darned much to pack in there, and events needed to evolve significantly over the last several releases, so the page grew and grew, so it's quite difficult to follow.

The first thing you need to do is to define a 'trigger', which tells when you want the event to occur.

In your case, the trigger will be when any of the derivatives go to zero, but there is a problem - the only variables that can be used to compute triggers are the variables that are output from dsolve (the working variables). So say you have the ODE:


When you call dsolve, the only variable that would appear for this problem is y(t):

> dsn := dsolve({D(y)(t)=sin(t),y(0)=-1},numeric):
> dsn(1);
                      [t = 1., y(t) = -0.540302207663555]

so the trigger can only depend on y(t) or t. The solution here, for detecting sign changes in y'(t), is to use the rhs, which is in terms of the data, so here you'd declare the trigger as follows:

> dsn := dsolve({D(y)(t)=sin(t),y(0)=-1},numeric,
>    events=[[sin(t),halt]]):
> dsn(1);
                      [t = 1., y(t) = -0.540302292330465]

> dsn(10);
Warning, cannot evaluate the solution further right of 3.1415926, event #1
triggered a halt
                [t = 3.14159265358979, y(t) = 1.00000000160744]

Does this help?

- Allan

This means that the discrete equations that are used in your dynamical system have a cyclic dependency.

For example, the equation that is used to define is given in terms of other variables including, and the equation that is used to define is given in terms of other variables including

Now discrete equations are solved using an iterative process, but the situation described by this warning indicates that we *may* be in a case where the iterative process will not converge, or may converge to an unexpected answer (i.e. changing the order of the defining equations may change the result of the iterative solution, as no single correct order could be determined).

You are correct though, the warning message should be improved, but it's difficult to describe this situation (and possible consequences) in sufficient detail in a warning message.

If you take your PDE2, and evaluate the right-hand-side to obtain the initial values for the time rate of change for S(R,t), you will see a curve that is extremely badly behaved around R=1.

Specifically it has a maximum around 7e15, and two minima around -3e-15.

The initial solution profile is a gaussian with a peak at R=1 of around 40.

What this means is that even with a really small time step on the order of 1e-14, part of the solution will become negative in the first time step.

This cannot be, as the PDE contains rational powers of S(R,t) (square roots, and powers of 3/2), so the solution rapidly becomes complex valued.

Also, the extremely rapid rate of change at the initial point is effectively a blow-up, so Float(infinity) rapidly enters the picture, and this is the source of the Float(undefined).

Perhaps a constant or a power is off in the formulation?


Allan Wittkopf

Maplesoft Math Developer



I have attached a worksheet (View on MapleNet or Download
View file details
) that shows one possible solution for the problem.

Basically, I identify the nonlinear equation and the nonlinear variable, then differentially extend the equation making it linear, also adding a constraint.

The problem can then be handled directly using the 'implicit=true' option.


Allan Wittkopf

Math Developer, Maplesoft

This is the same behavior that has been very recently observed in MaplePrimes for the rkf45 method.


Basically ODE's containing a solution component: x'(t)=x(t)/t which can be shiffted or added to another solution component, will be a problem for explicit RK techniques, as the solution is exact for any moderate or better order RK technique (see the excellent blog of Robert Israel on this issue in


The reason this is a problem is that a pair mechanism is used ot estimate error, and if both methods in the pair have zero error, the error estimate is always zero, which allows for an arbitrary step size.


Now combine this with the fact that the ODE form has a 1/t-type singularity, and that dsolve uses interpolation over the region of the step, and these odd results appear.


For ODE of this form, the other methods (rosenbrock, taylorseries, gear, etc.) work fine, so until this proiblem is resolved (and we are working on it) these methods should be used for problems of this class.


- Allan

The problem appears to be fixed in Maple 12, so I suspect you are using an older version of Maple?

I think this was a bug where the warning in the dsolve call was supressed, but not when using odeplot, and that was fixed during the development of Maple 12.

There is a separate issue with your example that appears in Maple 12 - specifically you form the ODEs with m,rho,CL as unassigned names, then assign them, and subsequently call dsolve from within a procedure (which does not evaluate the names). The problem is that then your ODE's appear to dsolve be parametric, and it fails in integrating them (as they appear to have additional unknowns m,rho,CL). You can remedy this by either:

1) Creating the ODE's after assigning m,rho,CL,etc.

2) Wrapping the ODE's in an eval in the call to dsolve as:



As a workaround for Maple 11 and the warnings, you could use:

_Env_in_maplet := true:

which is the mechanism that Maple uses to turn off warnings when dsolve/numeric and odeplot are being run from the interactive dsolve maplet.

- Allan

OK, now I have all the information.


The real issue is that the solution is non-polynomial, but both the 5th and 4th order solutions compute the answer exactly (thanks for the clue Robert!)


The method, as it is currently programmed, computes the solution from x0 to x0+h, then obtains an interpolant that should approximate the solution over x0..x0+h to 5th order accuracy. The problem is that since we get an exact solution for a non-polynomial input the interpolant is quite wrong, and the solution becomes incorrect for a region.


The really strange thing is that if you integrate past the one bad step, the solution becomes correct again (e.g. run the revised problem with initstep=1e-8, over t=0..10, and you will see a correct solution for 0..0.4 and 3.7...) which shows that it is purely a case where the interpolant is bad, but the solution points are being computed correctly.


- Allan

OK, this is a very interesting problem...

I have the exact coefficients for the RK pair, and for an equation as simple as x'=x/(1-t) I have obtained an exact form for the single step solution beginning at x(t0)=x0 for step size h as:


and the error estimate is indeed exactly zero.

The pair in use is the well known Fehlberg pair from

Klassische Runge-Kutta-Formeln vierter und niedrigerer Ordnung mit Schrittweiten-Kontrolle und ihre Anwendung auf Waermeleitungsprobleme, Computing, 6 (1970) 61-71

In contrast, the form of the error estimate for x'=-x is given by:

So here the error grows monotonically with the stepsize.

So this appears to be an issue with the method, and not the implementation.


- Allan Wittkopf

I have looked at this problem in detail, and the same issues can be observed for the simpler problem:

{diff(x(t),t)=x(t)/(1-t), x(0)=1}

which exhibits the same behavior for a singularity at t=1.


First of all, only the rkf45 method appears to have any difficulty with this problem, as choosing method=<anything else> appears to catch the singularity (or at least behave as though the solution is blowing up).


Tracing through the method, the problem is really that the step size chosen for the rkf45 method is dramatically increasing at each time step. The cause appears to be that the Runge-Kutta pair in use for this method dramatically underestimates the error - to such an extent that I believe that errors in equations of the form x'=x/(1-t) are in the effective nullspace of the pair in use.


Recall that the way that an rkf45 method works is to compute a pair of solutions of 4th and 5th order accuracy, then use their difference to estimate the error. The problem occurs when there is little or no difference between these solutions, but there is still significant error. I believe that this is one of those cases.


You can see that choosing a very small initial step size via 'initstep' causes the solution to be accurately computed over a (slightly) larger range.


This suggests the addition of a 'maxstep' option to restrict the maximum step size for a problem to avoid this difficulty.


- Allan Wittkopf

The default dsolve/numeric solvers (rkf45 and rosenbrock) have associated interpolants that can be used to compute the solution values. In addition, the procedure form of these solutions anticipates that additional solution points will be requested past the current solution point. With these two considerations in mind, the step size control in the schemes is such that the 'natural' step size, as dictated by the error control mechanism, is used in the integration. This means that when you request a solution at, say, t=2.0, the last step taken may be past t=2.0, and the t=2.0 value would be computed using the interpolant. This means that if you run the solver to obtain a value at t=10.0, you will always get the same value whether you request other intermediate time steps (e.g. 1.0, 2.0, ...) or not. This behavior can be modified by use of the 'steppast' option documented in ?dsolve,numeric,rkf45 as follows: > dprc := proc(n,t,y,yp) global last; > last := t; > yp[1] := y[2]; yp[2] := -y[1]; > end proc: > # Solution not using steppast > dsn1 := dsolve(numeric,method=rkf45,procedure=dprc, > start=0,initial=Array(1..2,[0,1]),procvars=[y,yp]); dsn1 := proc(x_rkf45) ... end proc > dsn1(2); [y = 2., y = 0.909297642748638890, yp = -0.416146922605180047] > last; 2.05387779892953137 > # Solution *using* steppast > dsn2 := dsolve(numeric,method=rkf45,procedure=dprc, > start=0,initial=Array(1..2,[0,1]),procvars=[y,yp], > steppast=false); dsn2 := proc(x_rkf45) ... end proc > dsn2(2); [y = 2., y = 0.909297635328363385, yp = -0.416146899596327335] > last; 2. - Allan
Implicit returns are a bit tricky, especially when you are dealing with an 'if' statement. The reason that your original piece of code did not work, is because the construct: ... friction:=-mu*N*(vn-k*v0)/(abs(vn)+k);#Placed here causes ERROR if vn=0 then friction:=RL_0 end if; end proc: will yield a 'NULL' as a return, as there is no clause in the 'if' that gives a value when vn<>0. Simpler example: > f := proc(x) > local friction; > > friction := x^2; > if x<1 then > friction := 1: > end if; > end proc: > > lprint(f(2)); NULL > lprint(f(0)); 1 Now the return of this procedure is being used in computing the function for the ODE, so if you get into a situation where vn<>0, the return will be NULL, and you will be attempting to do a computation with a NULL. This is easily fixed by putting in an explicit return, or as you observed, by adding in the second clause of the if, so that it always executes, and the relevant procedure always has a return: ... friction:=-mu*N*(vn-k*v0)/(abs(vn)+k);#Placed here causes ERROR if vn=0 then friction:=RL_0 end if; friction; end proc: or: ... friction:=-mu*N*(vn-k*v0)/(abs(vn)+k);#Placed here causes ERROR if vn=0 then friction:=RL_0 end if; return friction; end proc: should work. - Allan
1 2 3 Page 1 of 3