mathodediscrete-mathematicsmaplediscretization

Solving numerically a second order ODE with Maple


I am trying to solve the following ODE using the program Maple:

enter image description here

This ODE can easily be solved analytically on the interval [0,1/2] and the solution is given by

enter image description here

Now, using Maple, I got

restart;
ODE := diff(y(x),x$2)*x*abs(ln(x))-2*diff(y(x), x)=0:
initialPosition := y(0) = 0;
initialVelocity := D(y)(1/2) = 1;
solution := dsolve({ODE, initialPosition, initialVelocity})

which gives me the same solution I found earlier. However, when I try to find the solution numerically: solution_num := dsolve({ODE, initialPosition, initialVelocity}, numeric), the program displays the following error

enter image description here

As far as I understand, the way Maple discretizes the ODE together with the initial condition at 0, y(0) = 0, creates an error as the program cannot handle the singularity from the log(0) appearing in the equation. They suggest to use another type of discretization to circumvent this singular endpoint: the Midpoint Method. As I am just starting with Maple and I'm not comfortable with ODEs discretization, I was wondering if there was a function implementing such an algorithm to solve my problem. Any help would be appreciated.


Solution

  • restart;
    ODE:=diff(y(x),x$2)*x*abs(ln(x))
         -2*diff(y(x),x)=0:
    initPos:=y(0)= 0:
    initVel:=D(y)(1/2)=1:
    
    sol:=dsolve({ODE,initPos,initVel})
           assuming x>0,x<1;
    
        y(x) = (-x/ln(x)-Ei(1,-ln(x)))*ln(2)^2
    

    I put this style,etc on the plot of the exact solution only so that it can be nicely overlaid with the plot of the numeric solution.

    P1:=plot(eval(y(x),sol),x=0 .. 0.5,
             adaptive=false,numpoints=11,
             color=blue,style=point,
             symbol=solidcircle,symbolsize=20):
    

    Now, the numeric solution.

    numsol:=dsolve({ODE,initPos,initVel},numeric,
                   method=bvp[midrich],
                   maxmesh=10^4,abserr=5e-4):
    
    P2:=plot(X->eval(y(x),numsol(X)),0 .. 0.5):
    

    Now let's see them together.

    plots:-display(P1,P2);
    

    enter image description here

    Another (somewhat superior) way to plot that result returned by dsolve(...,numeric).

    plots:-odeplot(numsol,[x,y(x)],x=0 .. 0.5);
    

    enter image description here