@dharr 

My using isolate instead of solve is mainly a matter of laziness: isolate(equation, v) returns an equality that can directly be used in eval,  as solve (eq, v) requires writing v=solve (eq, v) to be used in eval.

@dharr 

I'm not at ease with eliminate and I often have to do several trials to get what I want.

Nice, I vote yup.

What is a(T) ?

@dharr 

Thanks

@shashi598

I adapted the Metropolis-Hastings' algorithm (developed to sample an unknown random variable) to find the values of alpha such that D(gamma)(xi__1) = 0.
(
More precisely, as the right boundary condition is gamma(xi__1)+xi__1* D(gamma)(xi__1) = 0, the previous condition is equivalent to gamma(xi__1) = 0, and it is equivalent to search for the values of alpha which make this identity true, or for values of alpha which minimize gamma(xi__1)^2.
I made this latter choice.
)

I got only 2 values (global minima of gamma(xi__1)^2) in the range alpha=1..10:

1. alpha = 3.7026096348364278322
2. alpha = 4.2985527068692070036

and 4 local minima (that is minima where gamma(xi__1)^2 is small but not null):

1. alpha = 5.3002873922249993805
2. alpha = 6.0544400684250000630
3. alpha = 7.5717744989250001680
4. alpha = 8.3303294592000002940

The type (local/global) of the minima are verified graphically.

You will note these values differ from those from the paper you cite.

Reply me if you still read this thread ad are interested to my code.

@Ronan 

Thank you Ronan

Could you check them?

@shashi598 

Explanation:

Is I'm not mistaken you want to find the value alpha such that diff(gamma(xi), xi) = 0 at xi=xi__1.
Which, give, that gamma(xi) + xi * diff(gamma(xi), xi) = 0 at xi=xi__1, is equivalent to find the value alpha such that gamma(xi) = 0 at xi=xi__1.

Here is the plot, with a logarithmic vertical axis, of gamma(xi__1)^2 in the range 1..6:


Then a focus on the range 3.5..4


and another on the range 4..6

Here are the values of alpha that minimize H(xi) = amma(xi__1)^2 when alpha is assumed to move in diffferent ranges:

This last plot completes the previous one by adding, in gray", the basins of attraction of red points, and in green those of the blue points.
Depending on your starting point in the range 4..6, NLPSolve will get anyone of the 5 points.
By default the starting point (if not specified with the option initialpoint) is the mid point of the searching range, thus 5 here.
Starting from 5 will return the first leftmost red point on this figure (that is your alpha2_result point).

Using DirectSearch may enable jumping out of this basin of attraction (the leftmost green rectangle) for it is not very deep.
Unfortunately the algorithm ends hitting the right boundary of the search domain (where the value of gamma(xi__1)^2 is lower than those at the red point, so it indeed minimizes the objective, at least in some way), but it misses the blue points because of their thinness and a bad tunning (I guess so) of the algorithm. 

Isn't it this a duplicate of thar previous question?

BTW, why do you persist in doing the same mistake (unknown=f(x) and parameter mu[f])?
Did you read my answer?


More of this your 4 equations are both identical except for the name of the dependent function, so their solutions will be fundamentaly identical it it is stupid to dsolve 4 equations.
The generic form of these equations is given by coeff(HPMeq, p, 1).
Before doing complex things  ask yourself if this equation can be solved formally; 

# example
test := coeff(HPMEq, p, 0);

# A solution is returned here
Without_bc := dsolve(test);   

# But not here
With_bc := dsolve({test, F[0](1) = 1, F[0](-1) = -1, (D(F[0]))(1) = 0, (D(F[0]))(-1)})

So trying to take the rhs of this last ouput generates the error message you got.

The only way is to dsolve numerically

{test, F[0](1) = 1, F[0](-1) = -1, (D(F[0]))(1) = 0, (D(F[0]))(-1)}

which implies to replace the many parameters by their numerical values.

Make an effort to understand all this and to submit us something that is worth us spending time answering, otherwise it's all a waste of  time for the people here.

@Saha 

 

Download eq1_only.mw

@shashi598 

Let's try to synchronize our exchanges.

What I propose is for you to read my previous reply (that I updated meanwhile), and for me to wait for your return.
Maybe, or maybe not, the procedure LaneEmdenSolution  my last reply contains will give you some answers.

If not, I propose you to explain again what you want to achieve . For instance:

• Is it enough for you to compute the solution (n=1) from 0 to Pi?
  If not, look below to the attached file.
   
• Is the problem for n <> 1 of the same nature (a condition, or boundary condition) at some specific value of xi?
  The LaneEmdenSolution procedure was written assuming this, but I can be cometely out of the way.

See you soon

BTW, here is a procedure that builds the numerical solution of the Lane-Emden equafion on any range 0..some_value:

>  restart Numerical solution of the Lane-Emden equation on arbitrary ranges Procedure:    1  Solve numerically the Emden equation  (result = EmdenN)    2  For the Lane-Emden equation         2.1   Use the known option option to instruct dsolve that it has to use the EmdenN solution         2.2   Solve numerically a boundary value problem for the Lane-Emden equation (result = LaneEmden_bc)                 from 0 to a given abscissa xi = BC_IC_transition_point with conditions        D(gamma)(0) = 0,        gamma(BC_IC_transition_point)+BC_IC_transition_point*D(gamma)(BC_IC_transition_point) = 0         2.3   Solve numerically an initial value problem for the Lane-Emden equation (result = LaneEmden_ic)                 from  xi = BC_IC_transition_point  to xi = EndPoint with conditions        gamma(BC_IC_transition_point)    = eval(gamma(xi), LaneEmden_bc(BC_IC_transition_point)        D(gamma)(BC_IC_transition_point) = eval(diff(gamma(xi), xi), LaneEmden_bc(BC_IC_transition_point)                  >  LaneEmden := proc(n, alpha, BC_IC_transition_point, EndPoint)   local gamma:   local ode1, ode2, EmdenN, LaneEmden_bc, LaneEmden_ic, f, bc, End_bc, ic:   ode1   := diff(xi^2*(diff(theta__0(xi), xi)), xi)/xi^2 = -theta__0(xi)^n;   EmdenN := dsolve({ode1, theta__0(0) = 1, (D(theta__0))(0) = 0}, theta__0(xi), numeric);   f := proc(z)     if not type(evalf(z),'numeric') then       'procname'(z);     else       eval(theta__0(xi), EmdenN(z));     end if;   end proc;   printf("Emden = 0 at xi = %a\n", fsolve(f(xi)=0, xi=0..EndPoint));   ode2         := diff(gamma(xi), xi, xi)-2*gamma(xi)/xi^2+alpha^2*gamma(xi) = -f(xi)^n*xi^2:   bc           := D(gamma)(0) = 0, gamma(BC_IC_transition_point)+BC_IC_transition_point*D(gamma)(BC_IC_transition_point) = 0;   LaneEmden_bc := dsolve({ode2, bc}, known=f, numeric, method=bvp[midrich], maxmesh=1024):   End_bc       := eval([gamma(xi), diff(gamma(xi), xi)], LaneEmden_bc(BC_IC_transition_point));   ic           := gamma(BC_IC_transition_point) = End_bc[1], D(gamma)(BC_IC_transition_point) = End_bc[2];   LaneEmden_ic := dsolve({ode2, ic}, known=f, numeric):      DocumentTools:-Tabulate(     [       plots:-display(         plots:-odeplot(           EmdenN           , [xi, theta__0(xi)]           , xi=0..EndPoint           , color=blue           , gridlines=true           , title = cat("Emden, n=", N, ", numerical solution")         )       ),              plots:-display(         plots:-odeplot(           LaneEmden_bc           , [xi, gamma(xi)]           , xi=0..BC_IC_transition_point           , color=blue           , gridlines=true           , title = cat("Lane-Emden, n=", N, ", numerical solution")         ),         plots:-odeplot(           LaneEmden_ic           , [xi, gamma(xi)]           , xi=BC_IC_transition_point..EndPoint           , color=red           , gridlines=true           , title = cat("Lane-Emden, n=", N, ", numerical solution")         )       ),              plots:-display(         plots:-odeplot(           LaneEmden_bc           , [xi, gamma(xi)+xi*D(gamma)(xi)]           , xi=0..BC_IC_transition_point           , color=blue           , gridlines=true           , title = typeset('gamma(xi)+xi*D(gamma)(xi)', `, n`=N)         ),         plots:-odeplot(           LaneEmden_ic           , [xi, gamma(xi)+xi*D(gamma)(xi)]           , xi=BC_IC_transition_point..EndPoint           , color=red           , gridlines=true           , title = typeset('gamma(xi)+xi*D(gamma)(xi)', `, n`=N)         )       )     ], width=80   ): end proc: >  LaneEmden(1, 1.4, Pi, 10) Emden = 0 at xi = 3.141592311   >  LaneEmden(3, 1.4, Pi, 10) Emden = 0 at xi = 6.896848025   >  LaneEmden(3, 1.4, 6.896848025, 10) Emden = 0 at xi = 6.896848025   >

Download General_Emden_and_Lane-Emden_2.mw

@shashi598

Here are the numerical solutions for the Emden and Lane-Emden equations when n=2.

I used the same condition at xi=Pi than for the case n=1.

Download The_case_n=2.mw

To complete this, you will find in this file General_Emden_and_Lane-Emden.mw a procedure named 
LaneEmdenSolution whose arguments are, in this order:

1. The vamlue of n.
2. The value of alpha.
3. The value of the right end point X where the right boundary condition 

   gamma(X) + X * D(gamma)(X) = 0

   is imposed

This procedure draws the numerical solutions of the Emden and Lane-Emden equations and the evolution of

gamma(xi)+xi*diff(gamma(xi), xi)  , in the range 0.. right_end_point

@shashi598  I'm sorry but I believe I deleted your last reply.

Fortunately I had another window opened with your text in it: could you copy/paste it in a new reply?

So after going through all your files which explains every step in detail. So an overview of the approach as i understand taken in your file:

- The Emden equation was solved numerically using dsolve to get EmdenN. 

- This was then used in the Lane-Emden equation ode2 to define the source term.

- An attempt was made to solve the Lane-Emden equation numerically using dsolve, but it failed due to a singularity at the initial point.

- To get around this, the initial point was shifted to ε > 0 using a parameter.

- However, there were still issues matching the numerical solution LaneEmdenN to the exact solution LaneEmdenE at ε and π.

- To improve the match, the initial conditions were tuned by defining a cost function Control that measures the error between LaneEmdenN and LaneEmdenE at π.

- Control was minimized over the initial conditions to find optimal values that provide a close match to LaneEmdenE.

- The minimized initial conditions were then used in dsolve to get an improved numerical solution LaneEmdenN_opt.

- LaneEmdenN_opt was shown to match LaneEmdenE to high precision, solving the original problem.

So in summary, the key steps were:

1) Shift initial point to avoid singularity 

2) Tune initial conditions to match known solution behavior

3) Minimize error between numerical and exact solutions

This allows dsolve to get an accurate numerical solution without needing the exact analytic form. The same general approach could be applied to other singular ODEs.

So, I want to know in general if I dont have access to Exact solution i.e EmdenE and as a result to LaneEmdenE which are used in defining ic2E to get the high accuracy, in cases for example n=3, whats the best strategy?. Should i evaluate at ApproximateXi__1 where Xi__1 is Pi when n=1?. But have limited accuracy since ic2E is not be available. Putting it in a simple way imagine a blackbox function that takes n ,alpha, epsilon and outputs D(gamma)(xi__1) .

And at last I think these tricks needs to be applied when we have singular ODEs. So, in case where its not the case a straight forward approach will work that doesn't require epsilons.right?

I really appreciate the time and effort you have put in your answer. I learned several new things from your file. Thanks again!

As a reply to the text above: you are right on all the points and understood perfectly what did.
nevertheless I believe the second approach detailed here SecondApproach is much simpler and suficient if you don't want to compute the numerical solution of the Lane-Emden equation beyond xi=Pi.

The "epsilon-eta-tau" strateby is tailored to your problem and I cannot guarantee you it will still work for other values of the exponent n.

If you have something specific in mind, for instance n <> 1, please let me know, I would be interested in working on it.

Again I'm sorry for having deleted your reply

@dharr 

Thank you dharr for the explanation.
I didn't know about workbook (it seems Maple 2015 doesn't have this heature, the only entry is  ExcelTools[WorkbookData] and  ExcelTools[WorkbookData] which refers to former).

Look these help pages

• Matrix
• LinearAlgebra
  more specificaly EigenValues and Eigenvectors