mmcdara

6214 Reputation

17 Badges

7 years, 315 days

MaplePrimes Activity


These are questions asked by mmcdara

There are things that seem simple but rapidly turn into a nightmare.

Here is an example: what I want is to the expression given at equation (4) in the attached file.

Using Int gives a wrong result.
Using int gives a right one but not of the desired form (some double integrals are nested while others are not).

I've been stuck on this problem for hours, can you please help me to fix it?

TIA

restart

use Statistics in
  # For more generality defina an abstract probability distribution.
  AbstractDistribution := proc(N)
    Distribution(
      PDF = (x -> varphi(seq(x[n], n=1..N)))
    )
  end proc:

  # Define two random variables pf AbstractDistribution type.
  X__1 := RandomVariable(AbstractDistribution(2)):
  X__2 := RandomVariable(AbstractDistribution(2)):

end use;

proc (N) Statistics:-Distribution(Statistics:-PDF = (proc (x) options operator, arrow; varphi(seq(x[n], n = 1 .. N)) end proc)) end proc

 

_R

 

_R0

(1)

F := (U1, U2) -> U1/(U1+U2);
T := mtaylor(F(X__1, X__2), [X__1=1, X__2=1], 2):

proc (U1, U2) options operator, arrow; U1/(U1+U2) end proc

(2)


Error: x[2] is droped out of the double integral in the rightmost term

use IntegrationTools in

J := eval([op(expand(T))], [seq(X__||i=x[i], i=1..2)]);
L := add(
       map(
         j ->  
         if j::numeric then
           j
         else
           (Expand@CollapseNested)(
             Int(
               j * Statistics:-PDF(X__1, x)
               , seq(x[i]=-infinity..+infinity, i=1..2)
             )
           )
         end if
         , J
       )  
     ):
ET := %
end use;

[1/2, (1/4)*x[1], -(1/4)*x[2]]

 

1/2+(1/4)*(Int(x[1]*varphi(x[1], x[2]), [x[1] = -infinity .. infinity, x[2] = -infinity .. infinity]))-(1/4)*x[2]*(Int(varphi(x[1], x[2]), [x[1] = -infinity .. infinity, x[2] = -infinity .. infinity]))

 

1/2+(1/4)*(Int(x[1]*varphi(x[1], x[2]), [x[1] = -infinity .. infinity, x[2] = -infinity .. infinity]))-(1/4)*x[2]*(Int(varphi(x[1], x[2]), [x[1] = -infinity .. infinity, x[2] = -infinity .. infinity]))

(3)


I want this

'ET' = 1/2
       +
       (1/4)*(Int(Int(x[1]*varphi(x[1], x[2]), x[1] = -infinity .. infinity), x[2] = -infinity .. infinity))
       -(1/4)*(Int(Int(x[2]*varphi(x[1], x[2]), x[1] = -infinity .. infinity), x[2] = -infinity .. infinity))

ET = 1/2+(1/4)*(Int(Int(x[1]*varphi(x[1], x[2]), x[1] = -infinity .. infinity), x[2] = -infinity .. infinity))-(1/4)*(Int(Int(x[2]*varphi(x[1], x[2]), x[1] = -infinity .. infinity), x[2] = -infinity .. infinity))

(4)


With int instead of Int one integral is double the other is double-nested

L := add(
       map(
         j ->  
         if j::numeric then
           j
         else
             int(
               j * Statistics:-PDF(X__1, x)
               , seq(x[i]=-infinity..+infinity, i=1..2)
             )
         end if
         , J
       )  
     ):
ET := %

1/2+int(int((1/4)*x[1]*varphi(x[1], x[2]), x[1] = -infinity .. infinity), x[2] = -infinity .. infinity)+int(-(1/4)*x[2]*(int(varphi(x[1], x[2]), x[1] = -infinity .. infinity)), x[2] = -infinity .. infinity)

(5)


As the expression of ET is now correct, I tried to use IntegrationTools to get the
form I want (equation (4)).

But as soon as I replace int by Int x[2] is again droped out.

So it's not even worth thinking about using CollapseNested!

 

use IntegrationTools in
  eval(ET, int=Int);  
end use;

1/2+Int(Int((1/4)*x[1]*varphi(x[1], x[2]), x[1] = -infinity .. infinity), x[2] = -infinity .. infinity)+Int(-(1/4)*x[2]*(Int(varphi(x[1], x[2]), x[1] = -infinity .. infinity)), x[2] = -infinity .. infinity)

(6)

 

Download Int_int.mw

A case where simplify(...) and simplify~(...) both return the wrong result.
Should we consider this a simplify bug?

restart:


A simple case

J := Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity)
     *
     Int(r[2]^2*varphi[2](r[2]), r[2] = -infinity .. infinity)

(Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity))*(Int(r[2]^2*varphi[2](r[2]), r[2] = -infinity .. infinity))

(1)

# OK

op(1, J) = simplify(op(1, J))

Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity) = Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity)

(2)

# OK

op(2, J) = simplify(op(2, J))

Int(r[2]^2*varphi[2](r[2]), r[2] = -infinity .. infinity) = Int(r[2]^2*varphi[2](r[2]), r[2] = -infinity .. infinity)

(3)

# But...
#
# Not OK

simplify(J)

(Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity))*(Int(r[1]^2*varphi[2](r[1]), r[1] = -infinity .. infinity))

(4)

# Not OK

simplify~(J)

(Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity))*(Int(r[1]^2*varphi[2](r[1]), r[1] = -infinity .. infinity))

(5)

# OK

map(simplify, J)

(Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity))*(Int(r[2]^2*varphi[2](r[2]), r[2] = -infinity .. infinity))

(6)


A slightly more complex case

J := (Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity))*(Int(r[2]^2*varphi[2](r[2]), r[2] = -infinity .. infinity))-(Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity))*(Int(r[2]*varphi[2](r[2]), r[2] = -infinity .. infinity))^2;

(Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity))*(Int(r[2]^2*varphi[2](r[2]), r[2] = -infinity .. infinity))-(Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity))*(Int(r[2]*varphi[2](r[2]), r[2] = -infinity .. infinity))^2

(7)

is(J=simplify(J))

false

(8)

is(J=simplify~(J))

false

(9)

is(J=map(simplify, J));
map(simplify, J);

false

 

(Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity))*(Int(r[1]^2*varphi[2](r[1]), r[1] = -infinity .. infinity))-(Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity))*(Int(r[1]*varphi[2](r[1]), r[1] = -infinity .. infinity))^2

(10)

add(map(u -> map(simplify, u), [op(J)]));

is(J=%);

(Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity))*(Int(r[2]^2*varphi[2](r[2]), r[2] = -infinity .. infinity))-(Int(r[1]^2*varphi[1](r[1]), r[1] = -infinity .. infinity))*(Int(r[2]*varphi[2](r[2]), r[2] = -infinity .. infinity))^2

 

true

(11)

 

Download Simplify_is_wrong.mw

Notional example:
I use mtaylor compute the Taylor expansion of a function f (U) of several variables U1, .., UN.
In the resul the terms are ordered this way:

  • the leftmost term is f (P) where P denotes the point where the expansion is done
  • followed by a succession of terms :
    • firstly ranked according to the total order of derivation of  f.
    • and among terms of same derivation order, ranked in some kind of lexicographic order

For instance

Vars    := [seq(U[i], i=1..2)]:
AtPoint := [seq(P[i], i=1..2)]:
mt      := mtaylor(f(Vars[]), Vars =~ AtPoint, 3)

 

f(P[1], P[2])+(D[1](f))(P[1], P[2])*(U[1]-P[1])+(D[2](f))(P[1], P[2])*(U[2]-P[2])+(1/2)*(D[1, 1](f))(P[1], P[2])*(U[1]-P[1])^2+(D[1, 2](f))(P[1], P[2])*(U[1]-P[1])*(U[2]-P[2])+(1/2)*(D[2, 2](f))(P[1], P[2])*(U[2]-P[2])^2

(1)

map(t -> op([0, 0], select(has, [op(t)], D)[]), [op(mt)][2..-1])

[D[1], D[2], D[1, 1], D[1, 2], D[2, 2]]

(2)
 

 

Download mtaylor.mw

How could I define another ordering of the terms in this mtaylor expansion?
For instance 1 being identified to some letter and 2 to another one such that 1 <  2 a lexicographic order would be 

D[1], D[1, 1], D[1, 2], D[2], D[2, 2]

Thanks in advance

Is it possible to transform relation (2) into relation (7) without using the hand-made relation (3) and the  Sum -> Int -> Sum trick?

restart

f := Product(x[i]^a*(1-x[i])^b, i)

Product(x[i]^a*(1-x[i])^b, i)

(1)

Lf := ln(f);

ln(Product(x[i]^a*(1-x[i])^b, i))

(2)

Sum(ln(x[i]^a*(1-x[i])^b), i)

Sum(ln(x[i]^a*(1-x[i])^b), i)

(3)

expand(%) assuming x[i] > 0, x[i] < 1, a > 0, b > 0

Sum(a*ln(x[i])+b*ln(1-x[i]), i)

(4)

eval(%, Sum=Int)

Int(a*ln(x[i])+b*ln(1-x[i]), i)

(5)

IntegrationTools:-Expand(%);

a*(Int(ln(x[i]), i))+b*(Int(ln(1-x[i]), i))

(6)

Lf = eval(%, Int=Sum)

ln(Product(x[i]^a*(1-x[i])^b, i)) = a*(Sum(ln(x[i]), i))+b*(Sum(ln(1-x[i]), i))

(7)

 

Download From2to7.mw

TIA

I recently answered a question concerning the Lane-Emden equation (see here LaneEmden) where the main topic was about finding its numerical solution.

The generic form of the Lane-Emden equation with parameter n is

LaneEmden := n -> (Diff(xi^2*(Diff(theta(xi), xi)), xi)) = -theta(xi)^n * xi^2

      d   /  2 / d            \\             n   2
n -> ---- |xi  |---- theta(xi)|| = -theta(xi)  xi 
      dxi \    \ dxi          //                  


I have just realized that I missed a "small" point in the original question: the OP ( @shashi598 ) wrote
"[...] Maple never comes out of evaluating [the] analytical solution when n=5 [...] ".
The important point here is that this solution (at least for some initial conditions) is known and simple (in the sense it doen't involve any special function).

So I tried for a few hours to verify this claim, and ended wondering myself if it might not be right?

Could you please tell me (I guess @shashi598 would be interested too in your return) if the differential equation LaneEmden(5) can be solved formally?
TIA.

Emden_equation.mw


EDITED:
After a little research it seems that very specigic method are used to build the analytic solution of the LaneEmden(n) (n not equal to 0, 1 and 5): serie expansions, homotopy, Adomian decomposition for instance.
I wasn't capable to find how the solution for LaneEmden(5) have been got for the first time (iseems to be atthe end of the 19th century).

2 3 4 5 6 7 8 Last Page 4 of 42