mmcdara

6635 Reputation

18 Badges

8 years, 129 days

MaplePrimes Activity


These are questions asked by mmcdara

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).

(I would prefer a solution for Maple 2015, but answers relative to newer versions are welcome)

Is there a simple way to force the result -y(1) + y(2) without using one of these two tricks?

# how can I get the expression of
int(diff(y(x), x), x=1..2);
                      / d                  \
                   int|--- y(x), x = 1 .. 2|
                      \ dx                 /

# Trick 1
int(diff(y(x), x), x);
eval(%, x=2)-eval(%, x=1)
                              y(x)
                          -y(1) + y(2)

# Trick 2
J := Int(diff(y(x), x), x = 1..2): 
value(IntegrationTools:-Parts(J, 1));
                          -y(1) + y(2)

TIA

4 5 6 7 8 9 10 Last Page 6 of 44