I found a bug when applying dsolve in the case when one of the ODE’s is deﬁned by means
of piecewise. Using piecewise is necessary when you need to provide initial conditions at \(x=0\)
when you have terms of the form `1/x`

(and you know there is no singularity). The following
instructions illustrate the problem:

> e1:=diff(mu(x),x)=x^2*exp(psi(x)); d 2 e1 := -- mu(x) = x exp(psi(x)) dx > e2:=diff(psi(x),x)=piecewise(x>0, -mu(x)/x^2, x=0,0); { mu(x) d { - ----- 0 < x e2 := -- psi(x) = { 2 dx { x { { 0 x = 0 > S:=dsolve({e1,e2,mu(0)=0,psi(0)=0},{mu(x),psi(x)}, type=numeric); Error, (in f) unable to store 'piecewise()' when datatype=float[8]

These are the equilibrium equations for an isothermal sphere (a standard text
book problem in newtonian astrophysics) and initial conditions should be given as
`mu(0)=0,psi(0)=0`

for physical reasons. However, since there is a term `mu(x)/x^2`

in the second equation, I normaly used piecewise to be able to set these initial
conditions.

This works very well in all previous Maple releases (from V5.1 to 7). I am using command line version of Maple8 for the Mac OSX but the same problem occurs in Windows and linux versions of Maple 8.

Is there a new feature of dsolve that allows you to set this type of initial conditions without using piecewise or is this a real bug?

It is corrected with Maple 9. (U. Klein)

Yes, this is indeed a bug.

There is, however, a workaround for this bug.

`dsolve/numeric`

allows a procedure-based speciﬁcation of a system. This means that you
can deﬁne the rhs of the ODE as a computation procedure to be used by `dsolve/numeric`

,
and allows much more ﬂexibility than just piecewise.

Information on this can be found in `?numeric,IVP`

.

The relevant options are `'procedure','procvars','initial', and 'start'`

.

For the same problem in Maple 8, we make the mapping:

Y[1] <-> mu(x), YP[1] <-> diff(mu(x),x) Y[1] <-> eta(x), YP[1] <-> diff(eta(x),x)

We can then deﬁne:

> pvars := [mu(x),psi(x)]; pvars := [mu(x), psi(x)] > dproc := proc(N,x,Y,YP) > YP[1] := x^2*exp(Y[2]); > if x>0 then > YP[2] := -Y[1]/x^2; > else > YP[2] := 0; > end if; > end proc; dproc := proc(N, x, Y, YP) YP[1] := x^2*exp(Y[2]); if 0 < x then YP[2] := -Y[1]/x^2 else YP[2] := 0 end if end proc > init := array([0,0]); init := [0, 0]

Where I note that I used an `'if'`

condition to avoid \(x=0\) for the equation for
`YP[2] (i.e. diff(eta(x),x))`

.

Now call `dsolve/numeric`

with this information:

> S := dsolve(numeric,procedure=dproc,procvars=pvars,initial=init,start=0); S := proc(x_rkf45) ... end proc

And things are identical to the prior versions:

Maple 8:

> S(1); [x = 1., mu(x) = 0.302901326325882625, psi(x) = -0.158827757643356910] > S(10); [x = 10., mu(x) = 25.1061282275108368, psi(x) = -3.73655959664163140]

Maple 7:

> e1:=diff(mu(x),x)=x^2*exp(psi(x)); d 2 e1 := -- mu(x) = x exp(psi(x)) dx > e2:=diff(psi(x),x)=piecewise(x>0, -mu(x)/x^2, x=0,0); { mu(x) d { - ----- 0 < x e2 := -- psi(x) = { 2 dx { x { { 0 x = 0 > S:=dsolve({e1,e2,mu(0)=0,psi(0)=0},{mu(x),psi(x)}, type=numeric); S := proc(rkf45_x) ... end proc > S(1); [x = 1., mu(x) = .302901367581909331, psi(x) = -.158827753673755290] > S(10); [x = 10., mu(x) = 25.1061294255031591, psi(x) = -3.73655980695852863]

Removing the piecewise part seems to do it:

> e1:=diff(mu(x),x)=x^2*exp(psi(x)): > e2:=diff(psi(x),x)=-mu(x)/x^2: > S:=dsolve({e1,e2,mu(0)=0,psi(0)=0},{mu(x),psi(x)}, type=numeric): > S(.3); [x = 0.3, mu(x) = 0.00891969113809897886, psi(x) = -0.0149329128389725352] > plots[odeplot](S,[[x,mu(x)],[x,psi(x)]],x=0..5);

Apart from the Maple problem, there’s a mathematical problem here. Putting a discontinuous term into a diﬀerential equation is likely to result in having no solutions, and you can’t expect Maple’s numerics to handle the delicate issues of analysis.

What I would do is ﬁnd a series solution near \(x=0\), evaluate it at some positive \(x\) (hopefully well within the radius of convergence), and use that as the starting point for a numerical solution.

In this case, after noting that according to the ﬁrst equation, in any solution that is analytic
at 0 we must have `mu(x) = O(x^3)`

:

> Order:= 17; e1:=diff(mu(x),x)=x^2*exp(psi(x)); e2:=diff(psi(x),x) = -mu(x)/x^2; S:= {mu(x) = add(a[i]*x^i, i=3..17), psi(x) = add(b[i]*x^i, i=1..15)}; se1:= series(eval(lhs(e1)-rhs(e1), S), x); se2:= series(eval(lhs(e2)-rhs(e2), S), x); eqs:= {coeffs(convert(se1,polynom),x), coeffs(convert(se2,polynom),x)}; solve(eqs); 2869 -629 168947 {b[12] = -----------, b[10] = ---------, a[15] = ------------, 13135122000 224532000 689593905000 61 -2869 629 b[8] = -------, a[13] = ----------, a[11] = --------, 1632960 1094593500 22453200 -61 -1 -1 a[9] = ------, a[7] = 1/315, a[5] = --, b[6] = ----, 204120 30 1890 -90151991 b[4] = 1/120, a[17] = ----------------, b[1] = 0, 3938960385360000 -168947 b[14] = -------------, a[4] = 0, b[3] = 0, b[2] = -1/6, 9654314670000 a[6] = 0, b[5] = 0, a[8] = 0, b[7] = 0, a[10] = 0, b[9] = 0, a[12] = 0, b[11] = 0, a[14] = 0, b[13] = 0, a[16] = 0, b[15] = 0, a[3] = 1/3} > Ssol:= subs(%, S); 2 4 6 61 8 Ssol := {psi(x) = -1/6 x + 1/120 x - 1/1890 x + ------- x 1632960 629 10 2869 12 168947 14 - --------- x + ----------- x - ------------- x , mu(x) 224532000 13135122000 9654314670000 3 5 7 61 9 629 11 = 1/3 x - 1/30 x + 1/315 x - ------ x + -------- x 204120 22453200 2869 13 168947 15 90151991 17 - ---------- x + ------------ x - ---------------- x } 1094593500 689593905000 3938960385360000

The radius of convergence seems to be well over 1, since the coeﬃcients are decreasing in
size. This series solution should be very accurate at, say,`x = 1/2`

.

> IC:= eval(Ssol,x=0.5); IC := {psi(0.5) = -0.04115395731, mu(0.5) = 0.04064923128} > Nsol:= dsolve({e1,e2, op(IC)},{psi(x),mu(x)}, numeric);

This numerical solution should be good for `x >= 1/2`

.

Thanks for your replies

Preben Alsholm’s remark (problem goes out by removing piecewise) is correct `**BUT ONLY**`

in Maple 8. If you try this system of ODE’s in previous releases you get error
messages

With Maple 7:

`"Error, (in S) cannot evaluate the solution further right of 0., probably a singularity"`

With Maple V5.1 and Maple 6:

`"Error, (in s) division by zero"`

That is why I solved this problem in earlier releases with piecewise and it worked very well.

Robert Israel’s remarks are absolutely relevant: in general you have to test the behavior of the involved terms near x=0 (or near any possible pole). However, in the case of the ODE system I mentioned (the isothermal sphere) I already knew that regularity at \(x=0\) is fulﬁled.

The key issue is that dsolve numerical should be able to handle piecewise functions as it did in previous releases. There should be (I guess) some Maple related problem in the example I mentioned. Am I reporting a Maple bug or (accidentaly) ﬁnding a Maple feature???