#### 7.52 bug in dsolve, Maple 8, and workaround (16.1.03)

##### 7.52.1 Roberto Sussman

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)

##### 7.52.2 Allan Wittkopf (30.1.03)

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]



##### 7.52.3 Preben Alsholm (30.1.03)

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);



##### 7.52.4 Robert Israel (30.1.03)

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;
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 .

##### 7.52.5 Roberto Sussman (3.2.03)

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???