[MUG] Bug in dsolve

Author: Maple User Group Posted: Thu, 13 Feb 2003 09:07:09 -0500

>> From: Maple User Group "maple_gr"

-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-

Date: Thu, 30 Jan 2003 15:24:28 -0800
From: Allan Wittkopf "awittkop"
Subject: bug in dsolve with piecewise
To: "maple-list"

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:
|\^/| Maple 7 (IBM INTEL LINUX)
._|\| |/|_. Copyright (c) 2001 by Waterloo Maple Inc.
\ MAPLE / All rights reserved. Maple is a registered trademark of
<____ ____> Waterloo Maple Inc.
| Type ? for help.
> 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);
bytes used=1029344, alloc=982860, time=0.08
[x = 10., mu(x) = 25.1061294255031591, psi(x) = -3.73655980695852863]

Hope this helps.

Sincerely,
Allan Wittkopf

-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-

Date: Fri, 31 Jan 2003 11:11:56 +0100
From: Preben Alsholm "P.K.Alsholm"
To: "maple-list" "sussman"
Subject: bug in dsolve with piecewise

--------------010103080602030207030205
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit

Removing the piecewise part seems to do it:

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

Preben Alsholm

sussman wrote:

>>>From: sussman "sussman"
>>>
>>>
>
>Hy MUG
>
>I found a bug when applying dsolve in the case when one of the ODE's is
>defined by means of piecewise. Using piecewsie 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:
>
> |\^/| Maple 8 (APPLE PPC OSX)
> ._|\| |/|_. Copyright (c) 2002 by Waterloo Maple Inc.
> \ MAPLE / All rights reserved. Maple is a registered trademark
>of
> <____ ____> Waterloo Maple Inc.
> | Type ? for help.
>
> > 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 equation 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?
>
>Roberto Sussman
>
>Dr Roberto A Sussman,
>Departamento de Gravitacion,
>Instituto de Ciencias Nucleares UNAM,
>Circuito Exterior CU,
>Mexico DF, 04510, MEXICO
>
>tels +52-555-6224690, 91, 92
>fax +52-555-6224693
>
>emails: "sussman"
> "sussky"
>
>web http://www.nuclecu.unam.mx/~sussman
>
>
>
>

--------------010103080602030207030205
Content-Type: text/html; charset=us-ascii
Content-Transfer-Encoding: 7bit

<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<title></title>
</head>
<body>
Removing the piecewise part seems to do it: 
 
> restart; 
> 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); 
 
Preben Alsholm 
 
sussman wrote: 
<blockquote type="cite"
cite= "mid200301302136.h0ULa2e19242"
<blockquote type="cite">
<blockquote type="cite">
<pre wrap="">From: sussman <a class="moz-txt-link-rfc2396E" href= "mailto:sussman"
</pre>
</blockquote>
</blockquote>
<pre wrap="">
Hy MUG

I found a bug when applying dsolve in the case when one of the ODE's is
defined by means of piecewise. Using piecewsie 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:

|\^/| Maple 8 (APPLE PPC OSX)
._|\| |/|_. Copyright (c) 2002 by Waterloo Maple Inc.
\ MAPLE / All rights reserved. Maple is a registered trademark
of
<____ ____> Waterloo Maple Inc.
| Type ? for help.

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

Roberto Sussman

Dr Roberto A Sussman,
Departamento de Gravitacion,
Instituto de Ciencias Nucleares UNAM,
Circuito Exterior CU,
Mexico DF, 04510, MEXICO

tels +52-555-6224690, 91, 92
fax +52-555-6224693

emails: <a class="moz-txt-link-abbreviated" href= "mailto:sussman"
<a class="moz-txt-link-abbreviated" href= "mailto:sussky"

web <a class="moz-txt-link-freetext" href="http://www.nuclecu.unam.mx/~sussman">http://www.nuclecu.unam.mx/~sussman</a>

</pre>
</blockquote>
 
</body>
</html>

--------------010103080602030207030205--

-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-

Date: Fri, 31 Jan 2003 10:30:14 -0800 (PST)
From: Robert Israel "israel"
To: "maple-list"
Subject: bug in dsolve with piecewise

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 .

Robert Israel "israel"
Department of Mathematics http://www.math.ubc.ca/~israel
University of British Columbia
Vancouver, BC, Canada V6T 1Z2

On Thu, 16 Jan 2003, sussman wrote:

> I found a bug when applying dsolve in the case when one of the ODE's is
> defined by means of piecewise. Using piecewsie 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:
>
> |\^/| Maple 8 (APPLE PPC OSX)
> ._|\| |/|_. Copyright (c) 2002 by Waterloo Maple Inc.
> \ MAPLE / All rights reserved. Maple is a registered trademark
> of
> <____ ____> Waterloo Maple Inc.
> | Type ? for help.
>
> > 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 equation 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?
>
> Roberto Sussman
>
> Dr Roberto A Sussman,
> Departamento de Gravitacion,
> Instituto de Ciencias Nucleares UNAM,
> Circuito Exterior CU,
> Mexico DF, 04510, MEXICO
>
> tels +52-555-6224690, 91, 92
> fax +52-555-6224693
>
> emails: "sussman"
> "sussky"
>
> web http://www.nuclecu.unam.mx/~sussman
>

-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-

Date: Mon, 3 Feb 2003 12:48:42 -0600
Subject: bug in dsolve with piecewise
To: Preben Alsholm "P.K.Alsholm"
From: roberto "sussky"

--Apple-Mail-2--581070289
Content-Transfer-Encoding: quoted-printable
Content-Type: text/plain;
charset=ISO-8859-1;
format=flowed

Thanks for your replies

Preben Alsholm's remark (problem goes out by removing piecewise) is=20
correct **BUT ONLY** in Maple 8. If you try this system of ODE's in=20
previous releases you get error messages

With Maple 7:
"Error, (in S) cannot evaluate the solution further right of 0.,=20
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=20
and it worked very well.

Robert Israel's remarks are absolutely relevant: in general you have to=20=

test the behavior of the involved terms near x=3D0 (or near any possible=20=

pole). However, in the case of the ODE system I mentioned (the=20
isothermal sphere) I already knew that regularity at x=3D0 is fulfiled.

The key issue is that dsolve numerical should be able to handle=20
piecewise functions as it did in previous releases. There should be (I=20=

guess) some Maple related problem in the example I mentioned. Am I=20
reporting a Maple bug or (accidentaly) finding a Maple feature???

Roberto Sussman

On Friday, January 31, 2003, at 04:11 AM, Preben Alsholm wrote:

> Removing the piecewise part seems to do it:
>
> > restart;
> > e1:=3Ddiff(mu(x),x)=3Dx^2*exp(psi(x)):
>
> > e2:=3Ddiff(psi(x),x)=3D-mu(x)/x^2:
> > S:=3Ddsolve({e1,e2,mu(0)=3D0,psi(0)=3D0},{mu(x),psi(x)}, =
type=3Dnumeric):
> > S(.3);
>
> =A0 [x =3D 0.3, mu(x) =3D 0.00891969113809897886,
>
> =A0=A0=A0=A0=A0=A0=A0 psi(x) =3D -0.0149329128389725352]
>
> > plots[odeplot](S,[[x,mu(x)],[x,psi(x)]],x=3D0..5);
>
> Preben Alsholm
>

On Friday, January 31, 2003, at 12:30 PM, Robert Israel wrote:

> From: Robert Israel "israel"
> Date: Fri Jan 31, 2003 12:30:14 PM America/Mexico_City
> To: "maple-list"
> Cc: "sussman"
> Subject: Re: [MUG] bug in dsolve with piecewise
>
>
> 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=3D0, 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) =3D O(x^3):
>
>> Order:=3D 17;
> e1:=3Ddiff(mu(x),x)=3Dx^2*exp(psi(x));
> e2:=3Ddiff(psi(x),x) =3D -mu(x)/x^2;
> S:=3D {mu(x) =3D add(a[i]*x^i, i=3D3..17), psi(x) =3D add(b[i]*x^i,=20=

> i=3D1..15)};
> se1:=3D series(eval(lhs(e1)-rhs(e1), S), x);
> se2:=3D series(eval(lhs(e2)-rhs(e2), S), x);
> eqs:=3D {coeffs(convert(se1,polynom),x),=20
> coeffs(convert(se2,polynom),x)};
> solve(eqs);
>
> 2869 -629 168947
> {b[12] =3D -----------, b[10] =3D ---------, a[15] =3D ------------,
> 13135122000 224532000 689593905000
>
> 61 -2869 629
> b[8] =3D -------, a[13] =3D ----------, a[11] =3D --------,
> 1632960 1094593500 22453200
>
> -61 -1 -1
> a[9] =3D ------, a[7] =3D 1/315, a[5] =3D --, b[6] =3D ----,
> 204120 30 1890
>
> -90151991
> b[4] =3D 1/120, a[17] =3D ----------------, b[1] =3D 0,
> 3938960385360000
>
> -168947
> b[14] =3D -------------, a[4] =3D 0, b[3] =3D 0, b[2] =3D =
-1/6,
> 9654314670000
>
> a[6] =3D 0, b[5] =3D 0, a[8] =3D 0, b[7] =3D 0, a[10] =3D 0, =
b[9] =3D 0,
>
> a[12] =3D 0, b[11] =3D 0, a[14] =3D 0, b[13] =3D 0, a[16] =3D =
0,
>
> b[15] =3D 0, a[3] =3D 1/3}
>
>> Ssol:=3D subs(%, S);
>
>
> 2 4 6 61 8
> Ssol :=3D {psi(x) =3D -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
> =3D 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=20
> coefficients
> are decreasing in size. This series solution should be very accurate=20=

> at,
> say, x =3D 1/2.
>
>> IC:=3D eval(Ssol,x=3D0.5);
>
> IC :=3D {psi(0.5) =3D -0.04115395731, mu(0.5) =3D 0.04064923128}
>
>> Nsol:=3D dsolve({e1,e2, op(IC)},{psi(x),mu(x)}, numeric);
>
> This numerical solution should be good for x >=3D 1/2 .
>
> Robert Israel "israel"
> Department of Mathematics http://www.math.ubc.ca/~israel
> University of British Columbia
> Vancouver, BC, Canada V6T 1Z2
>
>
>

My original message was

> sussman wrote:
>
> From: sussman "sussman"
>
>
>
> Hy MUG
>
> I found a bug when applying dsolve in the case when one of the ODE's =
is
> defined by means of piecewise. Using piecewsie is necessary when you
> need to provide initial conditions at x=3D0 when you have terms of the
> form 1/x (and you know there is no singularity). The following
> instructions illustrate the problem:
>
> |\^/| Maple 8 (APPLE PPC OSX)
> ._|\| |/|_. Copyright (c) 2002 by Waterloo Maple Inc.
> \ MAPLE / All rights reserved. Maple is a registered =
trademark
> of
> <____ ____> Waterloo Maple Inc.
> | Type ? for help.
>
> > e1:=3Ddiff(mu(x),x)=3Dx^2*exp(psi(x));
> d 2
> e1 :=3D -- mu(x) =3D x exp(psi(x))
> dx
>
> > e2:=3Ddiff(psi(x),x)=3Dpiecewise(x>0, -mu(x)/x^2, x=3D0,0);
> { mu(x)
> d { - ----- 0 < x
> e2 :=3D -- psi(x) =3D { 2
> dx { x
> {
> { 0 x =3D =
0
>
> > S:=3Ddsolve({e1,e2,mu(0)=3D0,psi(0)=3D0},{mu(x),psi(x)}, =
type=3Dnumeric);
> Error, (in f) unable to store 'piecewise()' when datatype=3Dfloat[8]
>
> These are the equilibrium equation for an isothermal sphere (a =
standard
> text book problem in newtonian astrophysics) and initial conditions
> should be given as mu(0)=3D0,psi(0)=3D0 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?
>
> Roberto Sussman
>
> Dr Roberto A Sussman,
> Departamento de Gravitacion,
> Instituto de Ciencias Nucleares UNAM,
> Circuito Exterior CU,
> Mexico DF, 04510, MEXICO
>
> tels +52-555-6224690, 91, 92
> fax +52-555-6224693
>
> emails: "sussman"
> "sussky"
>
> web http://www.nuclecu.unam.mx/~sussman
>
>
>
>

--Apple-Mail-2--581070289
Content-Transfer-Encoding: quoted-printable
Content-Type: text/enriched;
charset=ISO-8859-1

Thanks for your replies=20

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=20

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. =20

Robert Israel's remarks are absolutely relevant: in general you have
to test the behavior of the involved terms near x=3D0 (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=3D0 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???

Roberto Sussman

On Friday, January 31, 2003, at 04:11 AM, Preben Alsholm wrote:

<excerpt>Removing the piecewise part seems to do it:

> restart;

> e1:=3Ddiff(mu(x),x)=3Dx^2*exp(psi(x)):

> e2:=3Ddiff(psi(x),x)=3D-mu(x)/x^2:

> S:=3Ddsolve({e1,e2,mu(0)=3D0,psi(0)=3D0},{mu(x),psi(x)}, =
type=3Dnumeric):

> S(.3);

=A0 [x =3D 0.3, mu(x) =3D 0.00891969113809897886,

=A0=A0=A0=A0=A0=A0=A0 psi(x) =3D -0.0149329128389725352]

> plots[odeplot](S,[[x,mu(x)],[x,psi(x)]],x=3D0..5);

Preben Alsholm

</excerpt>

On Friday, January 31, 2003, at 12:30 PM, Robert Israel wrote:

<excerpt>From: Robert Israel < "israel"

Date: Fri Jan 31, 2003 12:30:14 PM America/Mexico_City

To: "maple-list"

Cc: "sussman"

Subject: Re: [MUG] bug in dsolve with piecewise

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=3D0, 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) =3D O(x^3):

<excerpt>Order:=3D 17;

</excerpt> e1:=3Ddiff(mu(x),x)=3Dx^2*exp(psi(x));

e2:=3Ddiff(psi(x),x) =3D -mu(x)/x^2;

S:=3D {mu(x) =3D add(a[i]*x^i, i=3D3..17), psi(x) =3D add(b[i]*x^i,
i=3D1..15)};

se1:=3D series(eval(lhs(e1)-rhs(e1), S), x);

se2:=3D series(eval(lhs(e2)-rhs(e2), S), x);

eqs:=3D {coeffs(convert(se1,polynom),x),
coeffs(convert(se2,polynom),x)};

solve(eqs);

2869 -629 168947

{b[12] =3D -----------, b[10] =3D ---------, a[15] =3D ------------,

13135122000 224532000 689593905000

61 -2869 629

b[8] =3D -------, a[13] =3D ----------, a[11] =3D --------,

1632960 1094593500 22453200

-61 -1 -1

a[9] =3D ------, a[7] =3D 1/315, a[5] =3D --, b[6] =3D ----,

204120 30 1890

-90151991

b[4] =3D 1/120, a[17] =3D ----------------, b[1] =3D 0,

3938960385360000

-168947

b[14] =3D -------------, a[4] =3D 0, b[3] =3D 0, b[2] =3D -1/6,

9654314670000

a[6] =3D 0, b[5] =3D 0, a[8] =3D 0, b[7] =3D 0, a[10] =3D 0, =
b[9] =3D 0,

a[12] =3D 0, b[11] =3D 0, a[14] =3D 0, b[13] =3D 0, a[16] =3D 0,

b[15] =3D 0, a[3] =3D 1/3}

<excerpt>Ssol:=3D subs(%, S);

</excerpt>

2 4 6 61 8

Ssol :=3D {psi(x) =3D -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

=3D 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 =3D 1/2.

<excerpt>IC:=3D eval(Ssol,x=3D0.5);

</excerpt>

IC :=3D {psi(0.5) =3D -0.04115395731, mu(0.5) =3D 0.04064923128}

<excerpt>Nsol:=3D dsolve({e1,e2, op(IC)},{psi(x),mu(x)}, numeric);

</excerpt>

This numerical solution should be good for x >=3D 1/2 .

Robert Israel "israel"

Department of Mathematics http://www.math.ubc.ca/~israel

University of British Columbia

Vancouver, BC, Canada V6T 1Z2

</excerpt>

My original message was

<excerpt>sussman wrote:

<fixed><bigger>From: sussman
=
<underline><color><param>1997,1997,FFFD</param>< "sussman"
/color></underline>

</bigger></fixed>

<fixed><bigger>

Hy MUG

I found a bug when applying dsolve in the case when one of the ODE's
is=20

defined by means of piecewise. Using piecewsie is necessary when you=20

need to provide initial conditions at x=3D0 when you have terms of the=20=

form 1/x (and you know there is no singularity). The following=20

instructions illustrate the problem:

|\^/| Maple 8 (APPLE PPC OSX)

._|\| |/|_. Copyright (c) 2002 by Waterloo Maple Inc.

\ MAPLE / All rights reserved. Maple is a registered
trademark=20

of

<<____ ____> Waterloo Maple Inc.

| Type ? for help.

> e1:=3Ddiff(mu(x),x)=3Dx^2*exp(psi(x));

d 2

e1 :=3D -- mu(x) =3D x exp(psi(x))

dx

> e2:=3Ddiff(psi(x),x)=3Dpiecewise(x>0, -mu(x)/x^2, x=3D0,0);

{ mu(x)

d { - ----- 0 << x

e2 :=3D -- psi(x) =3D { 2

dx { x

{

{ 0 x =3D =
0

> S:=3Ddsolve({e1,e2,mu(0)=3D0,psi(0)=3D0},{mu(x),psi(x)}, =
type=3Dnumeric);

Error, (in f) unable to store 'piecewise()' when datatype=3Dfloat[8]

These are the equilibrium equation for an isothermal sphere (a
standard=20

text book problem in newtonian astrophysics) and initial conditions=20

should be given as mu(0)=3D0,psi(0)=3D0 for physical reasons. However,=20=

since there is a term mu(x)/x^2 in the second equation, I normaly used=20=

piecewise to be able to set these initial conditions. This works very=20

well in all previous Maple releases (from V5.1 to 7). I am using=20

command line version of Maple8 for the Mac OSX but the same problem=20

occurs in Windows and linux versions of Maple 8.

Is there a new feature of dsolve that allows you to set this type of=20

initial conditions without using piecewise or is this a real bug?

Roberto Sussman

Dr Roberto A Sussman,

Departamento de Gravitacion,

Instituto de Ciencias Nucleares UNAM,

Circuito Exterior CU,

Mexico DF, 04510, MEXICO

tels +52-555-6224690, 91, 92

fax +52-555-6224693

emails:
=
<underline><color><param>1997,1997,FFFD</param></co=
lor></underline>

=20
=
<underline><color><param>1997,1997,FFFD</param></color><=
/underline>

web
=
<underline><color><param>1997,1997,FFFD</param>http://www.nuclecu.unam.mx/=
~sussman</color></underline>

</bigger></fixed>

</excerpt>=

--Apple-Mail-2--581070289--

Previous by date: [MUG] Re: Problem with SWP Implementation of Maple, Robert Israel
Next by date: [MUG] Displaying maples graphs in latex., Theo H S Boafo
Previous thread: [MUG] What has MAPLE 7 against Mathematica 3?, Donohue
Next thread: [MUG] Displaying maples graphs in latex., Theo H S Boafo

Ready to buy?

Maple - single user licence

Add to shopping basket

$ 1,895.00

Upgrade to Maple 12 from v11

Add to shopping basket

$ 995.00

Upgrade to Maple 12 from v10 & below

Add to shopping basket

$ 1,395.00

Featured Downloads

What's New in Maple 11 for Professionals

Maple White Paper: Technical Knowledge - An Asset You Can Afford to Lose?

Maple in Electronics Application Pack

Maple in Robotics & Aerospace Application Pack

Maple in Finance Application Pack

Product Reviews

"Without the Maple software, we would have to spend weeks generating the equations of motion for every experiment. Then the chances that we did it right would basically be near zero. There would always be a mistake somewhere. It is very difficult to set up a dynamic motion model by hand."
- Jean-Claude PiedBeouf, Ph.D Manager of Robotics, Canadian Space Agency

"Its very good - highly accurate and easy to use. The speed of Maple allows me to change equations and quickly reintegrate them into the application, so more possibilities can be explored to achieve the precise effect desired."
Shawn Neely, Senior R & D Director for PDI/Dreamworks