I found a bug when applying dsolve in the case when one of the ODE’s is defined 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 specification of a system. This means that you
can define the rhs of the ODE as a computation procedure to be used by dsolve/numeric,
and allows much more flexibility 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 define:
> 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 differential 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 find 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 first 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 coefficients 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 fulfiled.
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) finding a Maple feature???