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:
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:
We can then deﬁne:
Where I note that I used an 'if' condition to avoid \(x=0\) for the equation for YP (i.e. diff(eta(x),x)).
Now call dsolve/numeric with this information:
And things are identical to the prior versions:
Removing the piecewise part seems to do it:
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):
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.
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???