Subroutine to solve nonlinear DE

Posted on 2009-06-24 01:51 By manjees (

47)

Categories:

How do I ...? with Maple

Hi all,

I have a non-linear differential equation

x := diff(y(r), r, r)+2*(diff(y(r), r))/r+ C*polylog(3/2, -exp(C2*(C3-y(r)))) = 0
with boundary values D(y)(0)=0 and y(infinity)=0.

The hint says to use a subroutine to solve polylog and then integrate the differential equation with help of a Runge-Kutta like scheme.

I was wondering if anyone could help me with this. I have no clue how to use subroutines and specially when y(r) is within polylog.

Thanks in advance,
MS

here C, C2 and C3 are just

Comment on 2009-06-24 01:56 from manjees (

47)

here C, C2 and C3 are just numbers.

Still at it...

Comment on 2009-06-24 02:44 from Robert Israel ( Maple Rating 4

2618)

You seem to be stuck on variations of this one problem. From what I told you in another thread, the only way to have y(infinity) = 0 in a solution of diff(y(r), r, r)+2*(diff(y(r), r))/r = f(y) with f continuous at 0 would be f(0) = 0. Now here you have f(y) = - C*polylog(3/2, -exp(C2*(C3-y))). But polylog(3/2, t) has no real zeros other than 0, so except in the trivial case C=0 your boundary condition is impossible.

to get you started

Comment on 2009-06-24 02:47 from PatrickT ( Maple Rating 4

558)

Help:

dsolve[numeric,bvp],
dsolve/numeric_bvp,advanced,

you've got the following problems: some undefined constants, the polylog function, a bvp defined by a second-order ode.

to get started, I tried this:

> C:=1: C2:=1: C3:=1:
> ode := diff(y(r), r, r) +2*(diff(y(r), r))/r + C*polylog(3/2, -exp(C2*(C3-y(r))));
> R := 100:
> bcs := D(y)(0)=0, y(R)=0:
> dsolve({ode, bcs}, {y(r)}, range=0..R, type=numeric):
Error, (in dsolve/numeric/bvp) system is singular at left endpoint, use midpoint method instead

don't know how to proceed from here, sorry.

you could also transform the bvp into a set of 2 odes, but I don't know if that would help.

Hello sir, here y(r) is a

Comment on 2009-06-24 02:55 from manjees (

47)

Hello sir,
here y(r) is a potential and y(r -> infinity)=0. this has a solution as our professor has solved this problem. this is a part of summer homework problem set.

Our Prof emailed a hint saying to use a subroutine to solve polylog and then integrate the differential equation with help of a Runge-Kutta like scheme.

Would you mind throwing some light on how to make a subroutine for this case.

Really appreciate your help
MS

@ Patrick

Comment on 2009-06-24 03:02 from manjees (

47)

Even if we transform the BVP into a set of two eqns, still y(r) is inside polylog ang it creates trouble.

Just in case you need proper values of constants, C=486, C1 = 1/43, C2=43....

I have tried all the basic known methods in maple. The correct way is to a use a subroutine to solve polylog and then integrate the differential equation with help of a Runge-Kutta like scheme (according to my prof).

Thank you for your time,
MS

taylor approximate polylog

Comment on 2009-06-24 06:10 from PatrickT ( Maple Rating 4

558)

If you're so confident that there's a solution, then why don't you start with an approximation of polylog in terms of polynomials. Check the definition of the polylog (in Maple's Help menu), take a few terms of the expansion and see if that gives anything.

Also note that if 3/2 is replaced by something else (say by 1), the polylog function sometimes simplifies to better known functions.

But if I were you I would listen very, very carefully to Robert's advice.

Regarding sol

Comment on 2009-06-24 06:31 from manjees (

47)

1) I have already tried substituting few polynomial terms in place of polylog but still the DE remains non-linear.

2) OK regarding Solution of the DE: if u approximate polylog by few polynomial terms in y(r), then substitute that polynomial in y(r)(polylog) as another variable say a(r), one can write the DE in terms of a(r) which is still non-linear. That nonlinear DE in a(r) can be solved by using 'equilibrium iteration method or equilibrium Newton's method' numerically. Then the solutions values can be interpolated to get an expression for a(r) and for y(r). 'Numerical methods for eng and scientists: by Hoffman'.

I have to use this value of y(r) in some other expressions and so i need accurate rather than approximate values. Thats why I am emphasizing on using a subroutine to evaluate polylog and then integrating the differential equation with the help of a Runge-Kutta like scheme (according to my prof).

Hope this helps and thank you
MS

nonlinear

Comment on 2009-06-24 07:32 from PatrickT ( Maple Rating 4

558)

It's highly unlikely that there exist transformations that can make this ode linear. Nonlinear it is, nonlinear it will certainly remain. You just don't have a choice, you have to solve it numerically, assuming a solution does exist.

>That nonlinear DE in a(r) can be solved by using 'equilibrium iteration method or equilibrium Newton's method' numerically.

why don't you show the Maple code you've used to do that?

Yeah i know, i cannot change

Comment on 2009-06-24 07:39 from manjees (

47)

Yeah i know, i cannot change it to a linear equation...

How can I post the maple worksheet with the code here......?

Subroutines and RK

Comment on 2009-06-24 07:53 from manjees (

47)

I will send you the worksheet Patrick whenever I know how to post it. I can also email it to you.

But still I am hanging in mid air.

In my prof's words regarding solving this nonlinear DE: "Solve without
approximations. They key ingredient is to use a subroutine to evaluate the
polylogs and then integrate the
differential equation with help of a Runge-Kutta like scheme. You can
find both in programs like Mathematica, Maple, Matlab etc. "

I will really appreciate if I can get some help to figure a way out.

Thanks
MS

Subroutines?

Comment on 2009-06-24 21:16 from Robert Israel ( Maple Rating 4

2618)

I have no idea why your prof thinks using a "subroutine" to evaluate the polylogs is the key ingredient.
Subroutines are what Fortran had instead of procedures. Maple already has the capability of evaluating
polylogs. A "Runge-Kutta like scheme" is used to solve an initial value problem. Again, Maple can do
this with dsolve(..., numeric), with rather better methods than Runge-Kutta. The real problem is that the
differential equation has no solutions satisfying the boundary condition y(infinity)=0. All solutions will
go to infinity as r -> infinity. I suggest telling that to your prof.

language

Comment on 2009-06-24 08:33 from PatrickT ( Maple Rating 4

558)

the meaning of "without approximations" seems to be improper here: numerical solutions of ode systems (such as Runge-Kutta schemes) are approximations. Did the prof actually write that or is it your recollection?

you can copy-paste the Maple code into the window. If it's long, try to break the execution blocks and to insert linebreaks to make it readable.

Before copy-pasting, remove the output (Edit->Remove Output) to remove images and the likes.

posting a worksheet is also very easy, once you've posted it you can just copy-paste the url link to the worksheet. Go into "my files" on the lhs menu.

By "without approx" my prof

Comment on 2009-06-24 08:56 from manjees (

47)

By "without approx" my prof meant not to approximate anything for polylog as i was taking few polynomial terms in place of polylog.

And here is the worksheet. You have to run this once and then delete this line ''' for i from 0 to no-1 do y[i+1] := (dr+i*dr)^2*(-2.731/20^2) end do''' from sheet and then u can set a for loop. This above line is for initial values. View

I guess, i need to evaluate the polylog function using evaluate command and insert this as a subroutine in Rk integration...If you have any clue about doing this, please let me know.

Again, thank you and appreciate your patience
MS

maple worksheet problem

Comment on 2009-06-24 21:12 from PatrickT ( Maple Rating 4

558)

This is the header of your worksheet. It's full of html formatting. Can you post a pure Maple version?

<!DOCTYPE html
PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">

@ Robert...about roots of polylog

Comment on 2009-06-24 21:38 from manjees (

47)

I see what you are saying sir. Its polylog(3/2, -exp(1/43*(43-y(r))).

Notice if polylog is given by Ln(a,z) then z is -exp(...) and not exp(..).

I guess thats why you are getting complex roots and I hope this makes things clear.

Again any advice is appreciated.

Thanks
MS

roots

Comment on 2009-06-24 22:57 from Robert Israel ( Maple Rating 4

2618)

No, it doesn't clarify anything. I was not missing the negative sign. The function polylog(3/2, z) is increasing on the interval [-1,1] and has no real roots other than z=0.

about BC's

Comment on 2009-06-25 15:47 from manjees (

47)

yeah, i asked my prof about it... He told me to first solve this DE at initial conditions y(0)=0 and D(y)(0)=0 .

Then calculate the asymptotic potential y(r->∞) and subtract this from the entire potential function.

So, the immediate problem reduces to soving this nonlinear DE at IC's y(0)=0 and y'(0)=0.

Thanks for this insight Robert. Any suggestions now on solving this DE for initial conditions.

Thanks again
MS

IC's

Comment on 2009-06-26 01:50 from Robert Israel ( Maple Rating 4

2618)

OK, now you may have a doable problem.

With the initial conditions y(0)=0, D(y)(0) = 0, even though the DE is singular there, Maple is able to find a series solution.

> de:= diff(y(r),r,r)+2*diff(y(r),r)/r+486*polylog(3/2,-exp(1-1/43*y(r))) = 0;
  ys:= dsolve({de, y(0) = 0, D(y)(0) = 0}, y(r), series, order=10);

> yn := convert(rhs(evalf(ys)),polynom);

= 127.6269029*r^2-74.07550960*r^4+43.19136809*r^6-22.65830305*r^8

This series appears to have radius of convergence at least 1. To play it safe, I'll use this for 0 < r < 0.1. At r = 0.1 we have:

> yn1:= y(0.1) = eval(yn, r = 0.1);
  dyn1:= D(y)(0.1) = eval(diff(yn,r), r = 0.1

For larger values of r, we may use these as initial values for a numerical solution.

> S:= dsolve({de, yn1, dyn1}, y(r), numeric);

> with(plots): odeplot(S, [r, y(r)], r = 0.01 .. 100);

For large r, this looks approximately logarithmic. In fact, I suspect that all solutions are asymptotic to 86*ln(r) + 43 - 43*ln(43/243) as r -> infinity

Robert's method also works with reals

Comment on 2009-06-26 01:57 from PatrickT ( Maple Rating 4

558)

you can replace Rober't 486 by your 9*(16*43)*Pi^2*sqrt(Pi/(2*(1/43)^3))/(16*Pi^2*sqrt(43)^3), which by the way simplifies to 387/2*2^(1/2)*Pi^(1/2) and it will work. The output is just messier that's all.

simplify(9*(16*43)*Pi^2*sqrt(Pi/(2*(1/43)^3))/(16*Pi^2*sqrt(43)^3));

division by zero

Comment on 2009-06-25 22:31 from PatrickT ( Maple Rating 4

558)

irrespective of whether there is or isn't a solution to your problem, note that you are dividing by zero:

diff(y(r), r, r) + 2*(diff(y(r), r))/r + ...

then you have your boundary conditions defined for r=0. Problem.

The following produces something. Not sure if it's what you're after, but it may be a starting point.

> restart:
   with(plots):
   with(DEtools):

> C:=486: C2:= 1/43: C3:=43:

> ode := diff(y(r), r, r) + 2*(diff(y(r), r))/r 
         + C*polylog(3/2, -exp(C2*(C3-y(r))));
   ini := D(y)(1)=1, y(1)=0:

> S := dsolve({ode, ini}, {y(r)}, range=1..10, type=numeric, 
       abserr=1e-25, output=listprocedure):

> p := plots[odeplot](S): plots[display]({p});

You can also set up the problem as follows:

> ode:= { diff(y(r), r) = z(r),
                diff(z(r), r) = -2*z(r)/r 
                                - C*polylog(3/2, -exp(C2*(C3-y(r)))) 
              };
> ini:= {z(1)=0, y(1)=0};

Note also that because of the presence of r, your system is non-autonomous, so you cannot do a phase diagram analysis. But you could divide by some fixed number such as 1 or 10000 and do a phase diagram analysis to see what's going on at a fixed value of r.

Thanks for this. I have two

Comment on 2009-06-25 22:41 from manjees (

47)

Thanks for this. I have two questions.

You have changed your initial conditions at 1 rather than at 0. I don't think we can do that as we don't know the value of D(y)(1).

I just tried this simple thing:

de := diff(y(r), r, r)+2*(diff(y(r), r))/r+9*(16*43)*Pi^2*sqrt(Pi/(2*(1/43)^3))*polylog(3/2, -exp(1/43*(43-y(r))))/(16*Pi^2*sqrt(43)^3) = 0; dsolve[interactive]({de, y(0) = 0, (D(y))(0) = 0})

and used truncated series method in 'solve symbolically' option.
I am getting a solution like this:

y(r) = -(3/4)*sqrt(86)*sqrt(Pi)*polylog(3/2, -exp(1))*sqrt(43)*r^2-(1161/80)*Pi*polylog(1/2, -exp(1))*polylog(3/2, -exp(1))*r^4+O(r^6)

I am not sure if that is correct though. Would you mind having a look.

Thanks again
MS

> I don't think we can do

Comment on 2009-06-26 01:33 from PatrickT ( Maple Rating 4

558)

> I don't think we can do that as we don't know the value of D(y)(1).

you may need to *find* that value by thinking about the original problem, but what seems clear from the information you give is that dividing by r while at the same time setting r=0 will cause problems. Perhaps you can change your variables such that w:= r-> diff(y(r), r))/r; and have the division by zero disappear, perhaps...

Does the following plot look like the solution you're looking for?

y:= r-> -(3/4)*sqrt(86)*sqrt(Pi)*polylog(3/2, -exp(1))*sqrt(43)*r^2
        -(1161/80)*Pi*polylog(1/2, -exp(1))*polylog(3/2, -exp(1))*r^4;
p:= plot(y(r),r=0..1.5): display(p);

Got plots

Comment on 2009-06-26 03:21 from manjees (

47)

Thanks a million. I have 2 questions now.

1) How can i get an expression for y(r) from the plot above?

2) How can i calculate asymptotic potential for the diff equation at r-> as i need to subtract it from the expression in step(1)?

Thanks again for your help.
MS

Generic interpolation

Comment on 2009-06-26 11:47 from manjees (

47)

I got y(r) values from the formula at few points n did an Interpolation to get a solution. But I am not able to write a generic code for doing this. So, i just copy pasted points (a stupid laborious process).

1) Is there any generic way to get an expression for y(r) from the plot above?

2) How can i calculate asymptotic potential for the diff equation at r->infinity as i need to subtract it from the potential expression in step(1) to get a complete potential that satisfies BC's ?

Thanks
MS

question unclear

Comment on 2009-06-26 18:21 from PatrickT ( Maple Rating 4

558)

> Is there any generic way to get an expression for y(r) from the plot above?

I don't understand -- show your code and it may become clearer.

The plot is constructed from:

y(r) = series((-3/4*86^(1/2)*Pi^(1/2)*polylog(3/2,-exp(1))*43^(1/2))*r^2+(-1161/80*Pi*polylog(1/2,-exp(1))*polylog(3/2,-exp(1)))*r^4+(-27/448*86^(1/2)*Pi^(3/2)*43^(1/2)*polylog(-1/2,-exp(1))*polylog(3/2,-exp(1))^2-81/2240*86^(1/2)*Pi^(3/2)*43^(1/2)*polylog(1/2,-exp(1))^2*polylog(3/2,-exp(1)))*r^6+(-15093/8960*Pi^2*polylog(-1/2,-exp(1))*polylog(3/2,-exp(1))^2*polylog(1/2,-exp(1))-387/512*Pi^2*polylog(-3/2,-exp(1))*polylog(3/2,-exp(1))^3-3483/17920*Pi^2*polylog(1/2,-exp(1))^3*polylog(3/2,-exp(1)))*r^8+O(r^10),r,10)

and:

y(r) = 127.3728494*r^2-73.78089457*r^4+42.93395171*r^6-22.47842706*r^8

you can plot transformations of y(r) by using the data stored in the solution S:

plots[odeplot](S, [r, 2*y(r)], r = 0.01 .. 100);

> How can i calculate asymptotic potential

I don't know what an "asymptotic potential" is, sorry.

Code

Comment on 2009-06-26 18:40 from manjees (

47)

First we solve DE for initial conditions and then get values of y(.1) and D(y)(.1). Then we solve DE for those conditions. It's basically combination of Robert's and your advice.

ade:= diff(y(r),r,r)+2*diff(y(r),r)/r+486*polylog(3/2,-exp(1-1/43*y(r))) = 0;
ys:= dsolve({de, y(0) = 0, D(y)(0) = 0}, y(r), series, order=10);

yn := convert(evalf(ys), polynom);yn1 := eval(yn, r = .1); dyn1 := eval(diff(yn, r));

dyn2 := eval(dyn1, r = .1);S := dsolve({de, y(.1) = 1.266393114, (D(y))(.1) = 25.18200436}, y(r), numeric);with(plots); odeplot(S, [r, y(r)], r = 0.1e-1 .. 200);
# I get a plot from this step, but I need an expression for y(r). A crude method is to do S(1);S(2) and so on and copy paste those values for interpolation.

Would you mind helping me with a generic code for that step? And asymptotic potential is the solution for the DE at r-> infinity.

Thank you
MS