## 5635 Reputation

7 years, 189 days

@AHSAN

## @jalal Lok to ImageTools;-Read help...

@jalal

First syntax:

Read( file, img, opts )


A few lines below it's written

file    -   string; the pathname of the image file to read


Ask yourself the question: is  this:///Images/Maple.jpg a string?
Obviously not, so Read will fail (for the moment your "Read" is cornered by the special symbol ":" , see the its help page).

Now, is "this:///Images/Maple.jpg" a valid path?
I assume you use Windows and that this: is a the name of some disk? Unfortunately Maple doesn't accept this ,otation: you must enter the full path of Maple.jpg, not a shortcut.
A full path is something like "//../../../Images/Maple.jpg", to get it:

• go to the directiory Images,
• copy paste the name full path which should appear in the address bar,
• paste it in your worksheet,
• enclose all this by double quotes,
• change all the backslashes by slashes (maybe it's no longer necessary in recent Maple versions?).

Now you have a correct full path.

## @jalal I thought you would have und...

I thought you would have understood that my image had full path

cat(currentdir(), "/Desktop/Maple.jpg"):

I f you have some image whose full path is My_Image_is_here, just do

img := Read(My_Image_is_here):

I updated my previous answer and attached the mw file.

## A much simpler, while limited, approach...

UPDATED

You will find in this new file another point of view to solve numerically the Emden and Lane-Emden equations.
This is simpler than what was done before but (is it a drawback for you?) the numerical solution of the Lane-Emden equation can be computed ONLY in the range 0..Pi

ANALYSIS
Looking to the exact solution (n=1) of the Lane-Emden equation shows it contains the integration constant _C2 when initial conditions

gamma(0) = 0, D(gamma)(0) = 0

are used.
Said otherwise these IC do not determine completely the solution, so the idea to write instead (epsilon is a priori a small strictly positive value):

gamma(epsilon) = 0, D(gamma)(epsilon) = 0

Unfortunately this leads to a  solution which doesn't verify

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

Thus the quite complex procedure I used to identify two extra parameters eta and tau in this new Initial Value Problem:

ic := gamma(epsilon) = eta, D(gamma)(epsilon) = tau;
sol := dsolve({ode2(1), ic}), ....):

# find eta and ray such that
sol(Pi) + Pi * D(sol)(Pi) = 0

CONSEQUENCE
The other approach presented here starts from the observation that

gamma(0) = 0, D(gamma)(0) = 0

does not define a unique solution of the Lane-Emden equation.
So the idea is now to replace these initial conditions by these boundary conditions (see the attached file)

bc := D(gamma)(0) = 0, gamma(Pi) + Pi * D(gamma)(Pi) = 0

The correctness of this approach can be verified by doing this

rhs(dsolve({ode2E(1), bc}));
(cos(xi alpha) xi alpha - sin(xi alpha)) Pi
-------------------------------------------
/     2    \                    2
xi \alpha  - 1/ sin(Pi alpha) alpha

/     /      2    \   2\
\-2 + \-alpha  + 1/ xi / sin(xi) + 2 xi cos(xi)
+ -----------------------------------------------
2            2
xi (alpha - 1)  (alpha + 1)


and by checking that this expression is equal to the one we get from the original approach (identification of _C2).

Once this understood, getting the numerical solution of the Lane-Emden equation (with a numerical solution of the Emden equation) becomes obvious (note that this require using the option method=bvp[midrich] in dsolve/numeric).
The whole procedure is much simpler than the one previously used but forbids computing the numerical solution of the Lane-Emden equation beyond xi=Pi.

Here is the file

 > restart

If you want to compute the numerical solution of the Lane-Emden equation ONLY in the range 0..Pi
you can do the things in a simpler way.

 > local gamma:

Exact solutions of Emden and Lane-Emden (n=1) equations.

 > ode1 := n -> (diff(xi^2*(diff(theta__0(xi), xi)), xi))/xi^2 = -theta__0(xi)^n;
 (1.1)
 > EmdenE := unapply( rhs(dsolve({ode1(1), theta__0(0) = 1, (D(theta__0))(0) = 0}, theta__0(xi))), xi);
 (1.2)
 > ode2E := n -> diff(gamma(xi), xi, xi)-2*gamma(xi)/xi^2+alpha^2*gamma(xi) = -EmdenE(xi)^n*xi^2: ode2E(1)
 (1.3)
 > LaneEmdenE := rhs(dsolve({ode2E(1), gamma(0) = 0, D(gamma)(0) = 0}));
 (1.4)
 > AtPi := eval(LaneEmdenE+xi*diff(LaneEmdenE, xi), xi=Pi);
 (1.5)
 > HereIs_C2 := isolate(AtPi, _C2);
 (1.6)
 > LaneEmdenE := eval(LaneEmdenE, HereIs_C2);
 (1.7)

The strategy I used in my previous file
(and the same of yours  if I understood correctly what you did)

 > LaneEmdenE := rhs(dsolve({ode2E(1), gamma(0) = 0, D(gamma)(0) = 0}));
 (1.1.1)
 > AtPi := eval(LaneEmdenE+xi*diff(LaneEmdenE, xi), xi=Pi);
 (1.1.2)
 > HereIs_C2 := isolate(AtPi, _C2);
 (1.1.3)
 > LaneEmdenE := eval(LaneEmdenE, HereIs_C2);
 (1.1.4)

A simpler strategy where the condition gamma(0) is discarded and
replaced by gamma(Pi)+Pi*D(gamma)(Pi) = 0.

(Note this substitution is going to transform the Initial Value Problem into a
Boundary Value Problem)

 > Alternative := rhs(dsolve({ode2E(1), D(gamma)(0) = 0, gamma(Pi)+Pi*D(gamma)(Pi) = 0}));
 (1.2.1)
 > simplify(LaneEmdenE-Alternative);
 (1.2.2)

What happens at xi=Pi?

 > LaneEmdenE := eval(LaneEmdenE, alpha=1.4);
 (1.3.1)
 > eval(diff(LaneEmdenE, xi), xi=Pi);
 (1.3.2)

Numerical solutions of Emden and Lane-Emden equations.

 > ode2 := n -> diff(gamma(xi), xi, xi)-2*gamma(xi)/xi^2+alpha^2*gamma(xi) = -f(xi)^n*xi^2;
 (2.1)
 > EmdenN := dsolve({ode1(1), theta__0(0) = 1, (D(theta__0))(0) = 0}, theta__0(xi), numeric);
 (2.2)
 > f := proc(z)   if not type(evalf(z),'numeric') then     'procname'(z);   else     eval(theta__0(xi), EmdenN(z));   end if; end proc;
 (2.3)
 > Alternative_ic := D(gamma)(0) = 0, gamma(Pi)+Pi*D(gamma)(Pi) = 0;
 (2.4)
 > LaneEmdenN := dsolve(eval({ode2(1), Alternative_ic}, alpha=1.4), known=f, numeric, method=bvp[midrich])
 (2.5)
 > DocumentTools:-Tabulate(   [     plots:-display(       plots:-odeplot(LaneEmdenN, [xi, gamma(xi)], xi=0..Pi, color=red, legend=numeric),       plot(eval(LaneEmdenE, alpha=1.4), xi=0..Pi, color=blue, legend=exact),       title = Lane-Emden     ),     plots:-odeplot(       LaneEmdenN, [xi, gamma(xi)+xi*diff(gamma(xi),xi)], xi=0..Pi, color=red,       legend=typeset(gamma(xi)+xi*diff(gamma(xi),xi)),       gridlines=true     )   ], width=60 )

Note that you cannot obtain the numecical solution of the Lane-Emden equation
right to xi=Pi:

 > plots:-odeplot(LaneEmdenN, [xi, gamma(xi)], xi=0..5, color=red, legend=numeric)
 >

 > # A boundary value problem solved formally rhs(dsolve({ode2E(1), D(gamma)(0) = 0, gamma(Pi)+Pi*D(gamma)(Pi) = 0}));
 (3.1)
 > simplify(%-Alternative);
 (3.2)

 > # The no-condition formal solution of the Lane-Emden solution
 > FreeLaneEmden := rhs(dsolve(ode2E(1)));
 (4.1)
 > # Serie expansions at xi=0 S0 := series(FreeLaneEmden , xi=0, 4); S1 := series(diff(FreeLaneEmden, xi), xi=0, 4);
 (4.2)
 > # FreeLaneEmden is finite at xi=0 iif _C1=0 LeftFixedLaneEmden := eval(FreeLaneEmden, _C1=0)
 (4.3)
 > # The series expansion of diff(LeftFixedLaneEmden, xi) at xi=0 # shows this quantity is equal to 0 WHATEVER the value of _C2. series(diff(LeftFixedLaneEmden, xi), xi=0, 2)
 (4.4)
 > # Thus _C1=0 implies that both FreeLaneEmden and its derivative will be # equal to 0 at xi=0. # So taking these initial conditions: ic := gamma(0) = 0, D(gamma)(0) = 0; # doesn't completely define the solution, which is why an integration constant # appears now: rhs(dsolve({ode2E(1), ic}));
 (4.5)

Which one of the two initial conditions at xi=0 must be replaced by
gamma(Pi)+Pi*D(gamma)(Pi)=0 to get a uniquely defined solution?

 > # Replace D(gamma)(0) # (no solution found) ic := gamma(0) = 0, gamma(Pi)+Pi*D(gamma)(Pi)=0; dsolve({ode2E(1), ic});
 (5.1)
 > # Replace gamma(0) # (no solution found) ic := D(gamma)(0) = 0, gamma(Pi)+Pi*D(gamma)(Pi)=0; dsolve({ode2E(1), ic});
 (5.2)
 >

## @shashi598 (I'm answering you f...

@shashi598

UPDATED

(I answer you from my personal account, not my professional one @sand15).

By the same procedure as before here is the file corresponding to the case alpha=1.4.
You will find also the explanation of the different values you get at xi=Pi (essentially a numerical precision issue, see the lines after the text  Check the value of f(xi)+xi*diff(f(xi), xi) at point xi=Pi) and a way to find the initial conditions at xi=epsilon such that the solution verifies

gamma(xi)+xi*diff(gamma(xi), xi)=0

at point  xi =Pi (see after the text Attempt to circumvent using the exact solution to define the initial conditions at xi=epsilon.).
The main idea is

• to relax the initial conditions by writting gamma(epsilon)=eta and D(gamma)(epsilon)=tau (the numerical solution of the Lane-Emden equation then is a function of xi parameterized by eta and tau)
I used
epsilon = 10-8

• and then to minimize some functional of eta and tau in order that
 gamma(xi ; eta, tau)+xi*diff(gamma(xi ; eta, tau), xi)=0  for xi=Pi

The optimal solution  gamma(xi ; etaopt, tauopt) is equal to the exact solution up to an error of the order of 10-6 (a value that could be improved if necessary).
I got the optimal values

eta = -2.410160197*10-17
tau = -3.174864904*10-11

You may observe their are extremely close to the original initial condition eta=tau=0 we would write for the limiting case  epsilon=0;  but these small variationse have significant consequences:

eta = tau = 0:
gamma(Pi) + Pi*D(gamma)(Pi) = -1.01987         (instead of 0)
gamma(Z) + Z*D(gamma)(Z)    = 0  if Z = 3.2267 (instead of Pi)

eta = -2.410160197*10-17 and tau = -3.174864904*10-11
gamma(Pi) + Pi*D(gamma)(Pi) = 3.223690001*10-14  

Here is the file (link to DropBox)  LaneEmden_1dot4.txt

LaneEmdenNumeric_1dot4.mw

@Carl Love

You wrote "I didn't expect my technique to be the fastest possible", yes I understood that.
My opinion is that your technique has the strong advantage to manage a list of columns of any longer, as my procedure would be quite complex to generalize to a list of length larger than 2 (a recursive formulation could be useful).

You are also right, your technique is "reasonably fast and intuitive", (It needs to compare the lebgth of our respective codes).

 > restart:
 > kernelopts(version)
 (1)
 > MultipleSC := proc(M)   local cols, col, SC, tri, e, c, r, s, k:   uses ListTools:      SC := (M, j) -> M[sort(M[.., j], output=permutation), ..];   cols := [_passed[2..-1]]:   if numelems(cols) > 2 then     error "sorting wrt more than 2 column not yet implemented"   end if:   tri := sort(M[.., cols[1]], output=permutation);   if numelems(cols) = 2 then     e := [{entries(M[.., cols[1]], nolist)}[]];     c := convert(M[.., cols[1]], list);     r := [            0,            PartialSums(              map(u -> Occurrences(u, c), e)            )[]          ];     s := NULL:     for k from 1 to numelems(e) do       s := s, tri[r[k]+1..r[k+1]][sort(T[tri[r[k]+1..r[k+1]], cols[2]], output=permutation)][];     end do;     tri := [s];   end if:   return M[tri]; end proc: # Car Love's procedure CL := proc(M, L) local r, c; (r,c):= op(1, M): M[sort[inplace](     rtable(1..r, i-> ArrayTools:-Alias(M, c*(i-1), [1..c])),     key= (R-> [seq](R[L])),     output= permutation )] end proc:
 > N := 10^5: T := <|>(   seq(     <|>(       LinearAlgebra:-RandomVector(N, generator=1..4)       , Vector(N, i -> StringTools:-Random(1, "ABCD"))     )     , i=1..20   ) ): CodeTools:-Usage( MultipleSC(T, 1, 2), iterations=10 ): CodeTools:-Usage( CL(T ,[1, 2]), iterations=10 ):
 memory used=57.10MiB, alloc change=184.39MiB, cpu time=480.20ms, real time=260.10ms, gc time=275.48ms memory used=115.26MiB, alloc change=320.00MiB, cpu time=1.97s, real time=1.12s, gc time=1.11s
 > N := 10^5: T := <|>(   seq(     <|>(       LinearAlgebra:-RandomVector(N, generator=1..4)       , Vector(N, i -> StringTools:-Random(2, 'alpha'))     )     , i=1..20   ) ): CodeTools:-Usage( MultipleSC(T, 1, 2), iterations=10 ): CodeTools:-Usage( CL(T ,[1, 2]), iterations=10 ):
 memory used=56.95MiB, alloc change=213.64MiB, cpu time=296.90ms, real time=199.80ms, gc time=74.06ms memory used=115.26MiB, alloc change=-152.60MiB, cpu time=2.07s, real time=1.20s, gc time=1.14s
 > N := 10^6: T := <|>(       LinearAlgebra:-RandomVector(N, generator=1..4)       , Vector(N, i -> StringTools:-Random(2, 'alpha')) ): CodeTools:-Usage( MultipleSC(T, 1, 2) ): CodeTools:-Usage( CL(T ,[1, 2]) ):
 memory used=280.21MiB, alloc change=95.49MiB, cpu time=2.26s, real time=1.48s, gc time=414.64ms memory used=0.84GiB, alloc change=124.71MiB, cpu time=21.11s, real time=11.91s, gc time=10.82s
 >

## For information, here are a few performa...

What follows in no way diminishes the interest of your solution (see the reply I sent you earlier).
I improved  my procedure MultipleSC to avoid building submatrices and and managing only indices.
Here are compared performances for MultipleSC and your proposal.

 > restart:
 > MultipleSC := proc(M)   local cols, col, SC, tri, e, c, r, s, k:   uses ListTools:      SC := (M, j) -> M[sort(M[.., j], output=permutation), ..];   cols := [_passed[2..-1]]:   if numelems(cols) > 2 then     error "sorting wrt more than 2 column not yet implemented"   end if:   tri := sort(M[.., cols[1]], output=permutation);   if numelems(cols) = 2 then     e := [{entries(M[.., cols[1]], nolist)}[]];     c := convert(M[.., cols[1]], list);     r := [            0,            PartialSums(              map(u -> Occurrences(u, c), e)            )[]          ];     s := NULL:     for k from 1 to numelems(e) do       s := s, tri[r[k]+1..r[k+1]][sort(T[tri[r[k]+1..r[k+1]], cols[2]], output=permutation)][];     end do;     tri := [s];   end if:   return M[tri]; end proc:
 > if false then T := <|>(   LinearAlgebra:-RandomVector(10, generator=1..4)   , Vector(10, i -> StringTools:-Random(1, "ABCD"))   , Vector(10, i -> convert(StringTools:-Random(1, "uvwxyz"), name)) ): T worted wrt col 1 = MultipleSC(T, 1), T sorted wrt col 1 and col 2 = MultipleSC(T, 1, 2), T sorted wrt col 2 and col 1 = MultipleSC(T, 2, 1): end if;
 > N := 10^5: T := <|>(   seq(     <|>(       LinearAlgebra:-RandomVector(N, generator=1..4)       , Vector(N, i -> StringTools:-Random(1, "ABCD"))     )     , i=1..20   ) ): CodeTools:-Usage( MultipleSC(T, 1, 2), iterations=10 ): CodeTools:-Usage( [seq](R[[1,2]])))[]> ):
 memory used=57.10MiB, alloc change=184.39MiB, cpu time=490.90ms, real time=267.00ms, gc time=276.44ms memory used=225.25MiB, alloc change=231.45MiB, cpu time=5.00s, real time=3.65s, gc time=1.84s
 > TL := convert(T, listlist): CodeTools:-Usage( [seq](R[[1,2]])))[]> ):
 memory used=187.76MiB, alloc change=0 bytes, cpu time=3.33s, real time=2.52s, gc time=1.08s
 > N := 10^5: T := <|>(   seq(     <|>(       LinearAlgebra:-RandomVector(N, generator=1..4)       , Vector(N, i -> StringTools:-Random(2, 'alpha'))     )     , i=1..20   ) ): CodeTools:-Usage( MultipleSC(T, 1, 2), iterations=10 ): TL := convert(T, listlist): CodeTools:-Usage( [seq](R[[1,2]])))[]> ):
 memory used=56.95MiB, alloc change=91.56MiB, cpu time=313.70ms, real time=213.60ms, gc time=78.80ms memory used=187.76MiB, alloc change=-90.06MiB, cpu time=3.22s, real time=2.48s, gc time=926.57ms
 > N := 10^6: T := <|>(       LinearAlgebra:-RandomVector(N, generator=1..4)       , Vector(N, i -> StringTools:-Random(2, 'alpha')) ): CodeTools:-Usage( MultipleSC(T, 1, 2) ): TL := convert(T, listlist): CodeTools:-Usage( [seq](R[[1,2]])))[]> ):
 memory used=280.21MiB, alloc change=71.07MiB, cpu time=2.26s, real time=1.44s, gc time=440.63ms memory used=1.76GiB, alloc change=105.98MiB, cpu time=30.01s, real time=22.78s, gc time=9.04s
 >

It is likely that one can reach even better performances with a more subtle coding.

## @Carl Love I've just seen your ...

You deserve a gold star.

## @sursumCorda Here is a way to sort ...

Here is a way to sort a matrix wrt to the order induced by the elements of some column c1, and next (if required) wrt to the order induced by the elements of an other column c2.
This is very partial work for

• the procedure MultipleSC is restricted to sort wrt 1 or 2 columns (a recursive formulation seems necessary for a larger number of sorting columns;
• the column arguments could be gathered into a list;
• there is no checking of the input types nor of the values of the colums;
• the sorting acts only on columns, not on rows;
• I didn't assess the efficienvy of this procedire;
• ...

On the partial hospital[1..5, 1..3] test case displayed here  sortrows/sortcols, MultipleSC seems to return the expected result.

 > restart:
 > MultipleSC := proc(M)   local cols, col, SC, S1, e, c, r, k:   uses LinearAlgebra, ListTools:      SC := (M, j) -> M[sort(M[.., j], output=permutation), ..];   cols := [_passed[2..-1]]:   if numelems(cols) > 2 then     error "sorting wrt more than 2 column not yet implemented"   end if:   S1 := SC(M, cols[1]):   if numelems(cols) = 2 then     e := [{entries(S1[.., cols[1]], nolist)}[]]:     c := convert(S1[.., cols[1]], list):     r := [            0,            PartialSums(              map(u -> Occurrences(u, c), e)            )[]          ]:    <,>(      seq(        SC(          SubMatrix(S1, [r[k]+1..r[k+1]], [1..-1])          , cols[2]        )        , k=1..numelems(e)      )    ):    S1 := %:   end if:   return S1; end proc:
 > T := <|>(   LinearAlgebra:-RandomVector(10, generator=1..4)   , Vector(10, i -> StringTools:-Random(1, "ABCD"))   , Vector(10, i -> convert(StringTools:-Random(1, "uvwxyz"), name)) ); T worted wrt col 1 = MultipleSC(T, 1), T sorted wrt col 1 and col 2 = MultipleSC(T, 1, 2), T sorted wrt col 2 and col 1 = MultipleSC(T, 2, 1)
 (1)
 >

REMARK: I always found this lack of column/row sorting was painful (with R, which I use to use, it is really easy to do this, likely as with Matlab).
I suggest that you submit a specific question on this subject.

## My mistake...

I realized as I woke up this morning "Time taken to sort 1000000 elements 300 times" and I was about to connect on Mapleprimes in order to remove my last claim as I read your reply.

There is no such things for matrices.
Nevertheless, as I use Maple 2015 right now, I advice you to give a look to DataFrames if tou use a more recent version: if something like sortrows/sortcols do exist it should be there.

For the record, to sort a matrix wrt to the order induced by some row or column I use to do this (but I guess you were capable to find this by yourself :-)

restart
# sort all the columns wrt the order induced by column 1
SC := (M, j) -> M[sort(M[.., j], output=permutation), ..];

# sort all the lines wrt the order induced by line 1
SR := (M,i) -> M[.., sort(M[i, ..], output=permutation)];

T := LinearAlgebra:-RandomMatrix(5, 3, generator=0..4):
SC(T, 1);
SR(T, 2);



## Performances...

@sursumCorda

X is a vector of 105 floating points, Y is the same vector converted into a list.
PartialSorting_perf.mw

Main results:

• sort@Trim and Trim both are almost identical (memory used/allocated, cpi/real/gc times.
• Operating on a vector is six times faster then operating on a list and uses three times less memory.
• The performances are almost independent on the number (K) of elements in the partial vectors/lists.
• The cpu times are extremely small compared to those obtained with Matlab

## @sursumCorda  Yeah, I noticed ...

Yeah, I noticed that.
What is strange is that if you execute

CodeTools:-Usage( sort(Trim(X, 0, 100*k/N)) )


before

CodeTools:-Usage( Trim(X, 0, 100*k/N) ):


you get reversed times and Trim is now logicallly faster than sort@Trim.
See Look_at_this.mw

## @Axel Vogt Yes, our answers interse...

@Axel Vogt

It is only after I submitted mine that I saw yours, I hope you don't mind?

## NO...

Integrating without numerical bounds mean means integrate formally and I doubt that a closed form expression does exist.
So you have to integrate numerically, which means there must be numerical bounds.

Look here: How_to_integrate_formal.mw