mmcdara

Here is a partial result...

Answered: mmcdara 7891

March 18 2024

3 1

See the edited answer for more details and explanations

>

restart:

>	with(Statistics):

>	A := [8, 9, 16, 18]: B := [78, 31, 59, 81]:

>	# Direct formula fior the mean Moyenne := add(A*~B) / add(B); evalf(%);

(1)

>	# Direct formula for the standard deviation # (1): biased estimator ET := sqrt(add((A -~ Moyenne)^~2*~B) / add(B)); evalf(%);

(2)

>	# Direct formula for the standard deviation # (2): unbiased estimator ET := sqrt(add((A -~ Moyenne)^~2 *~ B) / (add(B)-1)); evalf(%);

(3)

>	# Alternative way to compute the mean and the standard deviation (unbiased estimator) C := [seq(A[i]$B[i], i=1..numelems(A))]: [Mean, StandardDeviation](C);

(4)

Download Ecart_type.mw

A thirteen steps solution......

Answered: mmcdara 7891

March 16 2024

1 1

@Mariusz Iwaniuk @oggsait

... whose steps ar guided by the result the OP wants to get.
So it is not a solution that Maple produces by itself but a solution we force Maple to proivide...

In particular, step (5), I don't know how to force Maple to convert cosh(X)^2-1 into sinh(X)^2 without using eval (as I did) or simplify with the rule i use in eval.
I also used several steps to transform sqrt(2)*sqrt(p[2])/sqrt(p[1]) into sqrt(2*p[2]/p[1]) because I wasn't able to do this more simply using ad hoc assumptions.

To @oggsait : I don't know if you will find some interest in that?

>

restart:

>	with(IntegrationTools):

>	f := 1/(usqrt(1+pu^2/q))

1/(u*(1+p*u^2/q)^(1/2))

(1)

>	F := Int(f, u);

Int(1/(u*(1+p*u^2/q)^(1/2)), u)

(2)

>

`(1)` = value(F);
`(2)` = isolate(rhs(%)=X, u);
`(3)` = allvalues(rhs(%));
`(4)` = simplify([rhs(%)]) assuming p > 0, q < 0, x > x[0];
`(5)` = eval(rhs(%), cosh = (z -> sqrt(1+sinh(z)^2)));
`(6)` = eval(rhs(%), sinh(X)=S);
`(7)` = simplify(rhs(%)) assuming S > 0;
`(8)` = eval(rhs(%), S=sinh(X));
`(9)` = eval(rhs(%), [p=P^2, q=Q^2]);
`(10)` = simplify(rhs(%)) assuming P > 0, Q < 0;
`(11)` = eval(rhs(%), Q=P*sqrt(q/p));
`(12)` = eval(rhs(%), sinh(X) = 1/sech(X));
`(13)` = eval(rhs(%), [X=x-x[0], q=``(-2*p[2]), p=p[1]]); # assuming p[2] < 0

`(1)` = -arctanh(1/(1+p*u^2/q)^(1/2))

`(3)` = (u = (-(tanh(X)^2*q-q)/(tanh(X)^2*p))^(1/2), u = -(-(tanh(X)^2*q-q)/(tanh(X)^2*p))^(1/2))

`(4)` = [u = (q/sinh(X)^2)^(1/2)/p^(1/2), u = -(q/sinh(X)^2)^(1/2)/p^(1/2)]

`(5)` = [u = (q/sinh(X)^2)^(1/2)/p^(1/2), u = -(q/sinh(X)^2)^(1/2)/p^(1/2)]

`(6)` = [u = (q/S^2)^(1/2)/p^(1/2), u = -(q/S^2)^(1/2)/p^(1/2)]

`(7)` = [u = q^(1/2)/(p^(1/2)*S), u = -q^(1/2)/(p^(1/2)*S)]

`(8)` = [u = q^(1/2)/(p^(1/2)*sinh(X)), u = -q^(1/2)/(p^(1/2)*sinh(X))]

`(9)` = [u = (Q^2)^(1/2)/((P^2)^(1/2)*sinh(X)), u = -(Q^2)^(1/2)/((P^2)^(1/2)*sinh(X))]

`(11)` = [u = -(q/p)^(1/2)/sinh(X), u = (q/p)^(1/2)/sinh(X)]

`(13)` = [u = -(``(-2*p[2])/p[1])^(1/2)*sech(x-x[0]), u = (``(-2*p[2])/p[1])^(1/2)*sech(x-x[0])]

(3)

>

Download A_thirteen_steps_solution.mw

A solution...

Answered: mmcdara 7891

March 15 2024

0 7

Even if it is possible to usre Runge-Kutta schems to solbe boundary value problems, this is not allowable in Maple.
So either you write your own RK algorithm or you trust Maple to use the scheme it founds most suited to.

solution.mw

Explanation...

Answered: mmcdara 7891

March 08 2024

1 1

The "ciulprit" is the term

sqrt(gamma^2*sigma__d^2*sigma__e^4 + 4*sigma__e^2 + 4*sigma__v^2)

in the numerator of A.
The term under the radical is positive and can be noted T^2.
Now what is sqrt(T^2)?

b := sqrt(a^2):
simplify(b) assuming a > 0
                               a
simplify(b) assuming a < 0
                               -a

This is why you need assumptions.

Depending on the sign of T the numerator of A is equivalent to gamma^3 or gamma^5 when
gamma --> +oo (while the denominator of A is equivalent to gamma^4).

A detailed answer is developed in the attached file limits_signum_explained.mw

If I may offer some advice...

Answered: mmcdara 7891

March 07 2024

1 4

To get less confusing notations a good practice is to assign capital letters to tandom variables.
So change delta by Delta, nu by Nu and so on. The relations will be clearer.

There is no need no write nu__1-nu__01 where nu_1 is a random variable and nu__01 its mean. The only advantage in doing this is to add an unnecesary complexity!

Lastly, why not just ask Maple what those covariances are?

>

restart

>	with(Statistics):

>

RVS := [
         Nu__1    = RandomVariable(Normal(0, sigma__nu)),
         Nu__2    = RandomVariable(Normal(0, sigma__nu)),
         Delta__1 = RandomVariable(Normal(0, sigma__d)),
         Delta__2 = RandomVariable(Normal(0, sigma__d)),
         Delta__3 = RandomVariable(Normal(0, sigma__d3))
       ]

(1)

>	X__1 := beta__1(Nu__1+Nu__2)+alpha__1Delta__1+alpha__2sDelta__2; X__2 := beta__2(Nu__1+Nu__2)+alpha__2Delta__2+alpha__1sDelta__1; X__3 := beta__3(Nu__1+Nu__2)+alpha__3Delta__3;

(2)

>	A := X__1(-lambda__1X__1-lambda__1Delta__1+Nu__1); B := X__2(-lambda__2X__2-lambda__2Delta__2+Nu__2); C := X__3(-lambda__3X__2-lambda__3*Delta__3+Nu__1+Nu__2);

(3)

>	`Cov(A, B)` := Covariance(eval(A, RVS), eval(B, RVS)):

>	`Cov(A, B)` := simplify(`Cov(A, B)`);

(4)

>	# and so on for Cov(A, C) and Cov(B, C)

Download No_problem.mw

gamma l!!!...

Answered: mmcdara 7891

March 06 2024

2 5

One problem comes from the use of gamma, which is a protected name (Euler constant, see help page).
So when you write gamma, Maple will not consider it as an unknown in the solve process (which is why your solutions never contain something of the form gamma < (or >) something.

Another problem comes from using the assumptions directly in solve.
If you do this (look to my file for the notations), you will get solutions which contain b < 0

solve(expr > 0) assuming hyp:

# one non licit solution is 

{E = E, G = 1/(3*b), a = a, s = s, b < 0, -(-s+3*b)*E < c}

So you have to "filter" the solutions solve returns.

(
NOTES:

I used SOL_gt := solve(expr > 0) assuming hyp but SOL_gt := solve(expr > 0) seems to give the same result.
SOL_gt := solve(expr > 0, useassumptions=true) assuming hyp returns a wrong result (Maple 2015).
SOL_gt := solve({expr > 0, hyp}) seems to run endlessly (Maple 2015).

)

Here is the solution:

>

restart

>	# gamma is the Euler constant gamma; evalf(gamma);

(1)

>	expr := (-(4((Es-2c)b+Esb))bG^2+(8b^2E+4Esb+(4(Es-3a(1/2)-(1/2)c))b)G-bE-Es+2a-E(b+s))/(2(2b*G-1)^2);

(1/2)*(-4*((E*s-2*c)*b+E*s*b)*b*G^2+(8*b^2*E+4*E*s*b+4*(E*s-(3/2)*a-(1/2)*c)*b)*G-b*E-E*s+2*a-E*(b+s))/(2*G*b-1)^2

(2)

>	hyp := a > 0, b > 0, c > 0, a > 2c, G > 0, 2b*G > a/c, s > 0, E >= 0

(3)

>	SOL_gt := solve(expr > 0) assuming hyp: print~([%]):

${E = E, c = c, s = s, G < (1/2)/b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), b < 0}$

${E = E, c = c, s = s, G < (1/3)/b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), b < 0, (1/2)/b < G}$

${E = E, c = c, s = s, b < 0, (1/3)/b < G, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}$

${E = E, c = c, s = s, 0 < b, G < (1/3)/b, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}$

${E = E, c = c, s = s, 0 < b, G < (1/2)/b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), (1/3)/b < G}$

${E = E, c = c, s = s, 0 < b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), (1/2)/b < G}$

(4)

>

# Eliminate solutions that do not verify somple hypothesis of the form
# "0 < X" ot "0 <= X"

simplehyp := select(has, {hyp}, 0);

good_gt := NULL:

for z in {SOL_gt} do
  map(x -> (coulditbe(x) assuming simplehyp[]), [z[]]);
  if not(member(false, %)) then
    good_gt := good_gt, z
  end if:
end do:

good_gt := {good_gt}:

numelems({SOL_gt}), numelems(good_gt);

print(cat('_'$140)):

print~(good_gt):

${E = E, c = c, s = s, 0 < b, G < (1/3)/b, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}$

${E = E, c = c, s = s, 0 < b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), (1/2)/b < G}$

${E = E, c = c, s = s, 0 < b, G < (1/2)/b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), (1/3)/b < G}$

(5)

>	# I couldn't find a way to check if these four solutions verify the initial hypothesis (hyp)

>	SOL_lt := solve(expr < 0) assuming hyp: print~([%]):

${E = E, c = c, s = s, G < (1/2)/b, b < 0, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}$

${E = E, c = c, s = s, G < (1/3)/b, b < 0, (1/2)/b < G, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}$

${E = E, c = c, s = s, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1), b < 0, (1/3)/b < G}$

${E = E, c = c, s = s, 0 < b, G < (1/3)/b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1)}$

${E = E, c = c, s = s, 0 < b, G < (1/2)/b, (1/3)/b < G, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}$

${E = E, c = c, s = s, 0 < b, (1/2)/b < G, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}$

(6)

>

# Eliminate solutions that do not verify somple hypothesis of the form
# "0 < X" ot "0 <= X"

good_lt := NULL:

for z in {SOL_lt} do
  map(x -> (coulditbe(x) assuming simplehyp[]), [z[]]);
  if not(member(false, %)) then
    good_lt := good_lt, z
  end if:
end do:

good_lt := {good_lt}:

print(cat('_'$140)):

print~(good_lt):

${E = E, c = c, s = s, 0 < b, G < (1/3)/b, a < -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1)}$

${E = E, c = c, s = s, 0 < b, (1/2)/b < G, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}$

${E = E, c = c, s = s, 0 < b, G < (1/2)/b, (1/3)/b < G, -(4*E*G^2*b^2*s-4*G^2*b^2*c-4*E*G*b^2-4*E*G*b*s+G*b*c+E*b+E*s)/(3*G*b-1) < a}$

(7)

>	# check that no goo_gt solution is a good_lt solution and reciprovcally good_gt union good_lt: numelems(%); good_gt intersect good_lt: numelems(%);

(8)

Download Inequalities_mmcdara_1.mw

Why did you delete your previous questio...

Answered: mmcdara 7891

March 05 2024

0 0

Is it too much to ask that you respect those who gave you the answer?

A way, likely among many...

Answered: mmcdara 7891

March 05 2024

1 2

Download A_way.mw

Confusion between variance and raw momen...

Answered: mmcdara 7891

March 04 2024

2 2

The first point is : "How do you define the mean and variance of g?
Look here reply_0.mw and read the last comment at the end of my answer.

Assuming that you want expressions for Ey and Vy which depend on the thetas, there is no need to define the latters as random variables before computing expressions Ey and Vy.
Finally your formula for Vy does not represents the variance of g but its 2nd order raw moment:

>

restart:

>	with(Statistics):

If you expect expressions of Ey and Vy to depend upon the thetas you must not define them before computing these expressions.
Instead of what you will get two numbers.

>	x[1] := RandomVariable(Normal(theta[1], theta[3]));

>	x[2] := RandomVariable(Normal(theta[2], theta[4]));

(1)

>	g := 1 - (x[1] - 1)^2/9 - (x[2] - 1)^3/16;

(2)

>	MyEy := CodeTools:-Usage( expand(Mean(g)) ); MyVy := CodeTools:-Usage( expand(Variance(g)) );

memory used=6.08MiB, alloc change=32.00MiB, cpu time=126.00ms, real time=1.26s, gc time=0ns

memory used=6.97MiB, alloc change=0 bytes, cpu time=78.00ms, real time=204.00ms, gc time=0ns

(9/64)*theta[4]^4+(2/81)*theta[3]^4+(9/64)*theta[2]^2*theta[4]^4-(9/32)*theta[2]*theta[4]^4+(4/81)*theta[1]^2*theta[3]^2-(8/81)*theta[1]*theta[3]^2+(4/81)*theta[3]^2+(15/256)*theta[4]^6+(9/256)*theta[4]^2*theta[2]^4-(9/64)*theta[4]^2*theta[2]^3+(27/128)*theta[4]^2*theta[2]^2-(9/64)*theta[4]^2*theta[2]+(9/256)*theta[4]^2

(3)

>

# Compare to what you want

Ey := -(1/9)*theta[1]^2 + (2/9)*theta[1] + 137/144 - (1/16)*theta[2]^3 + (3/16)* theta[2]^2 - (3/16)*theta[2] - (3/16*(-1+theta[2]))*theta[4]^2 - (1/9)*theta[3 ]^2:

Vy := (1/20736)*(9*theta[2]^3 + 16*theta[1]^2 - 27*theta[2]^2 - 32*theta[1] + 27*theta[2] - 137)^2 + (15/256)*theta[4]^6 + (45/256*(theta[2]^2 - 2*theta[2] + 1))*theta[4]^4 + ( 1/768*(45*theta[2]^4 + 32*theta[1]^2*theta[2] - 180*theta[2]^3 - 32*theta[1]^2 - 64*theta[ 1]*theta[2] + 270*theta[2]^2 + 64*theta[1] - 436*theta[2] + 301))*theta[4]^2 + (1/27)*theta[3]^4 + (1/216)*theta[3]^2*(3*theta[2]^3 + 16*theta[1]^2 - 9*theta[2]^2 - 32*theta[ 1] + 9*theta[2] - 35) + (1/24)*theta[3]^2*(-1+theta[2])*theta[4]^2;

(4)

>	simplify(MyEy-Ey); simplify(MyVy-Vy);

(5)

You probably confounded the expression of the Variance with those of the Raw Moment of order 2:

>	M2 := Moment(g, 2);

(6)

>	simplify(M2-Vy);

(7)

>	VarianceFromM2 := M2 - MyEy^2: simplify(VarianceFromM2 - MyVy)

(8)

>

Download reply_1.mw

Last comment: neither x[1] nor x[2] are random variables (in a mathematical sense) if the parameters of teir distributions are random variables.
But, for instance, x[1] conditionnally to theta[1]=t[1] and theta[3]=t[3] is a gaussian random variables with parameters t[1] and t[3].
So do you want to compute conditional mean and variance or "whole" mean and variance?

If I am not mistaken...

Answered: mmcdara 7891

March 03 2024

1 2

A very likely non optimal code that could probably be enhanced.

>

restart

>

GeneratorsOf := proc(n)
  local N, F, G, s, t, r, d:

  uses combinat:

  N := 1:
  while fibonacci(N) < n do N := N+1 end do: N := max(6, N-1):

  F := {seq(fibonacci(i), i=2..N)};
  G := convert(convert~(powerset(F) minus {{}}, list), list):
  s := select((g -> add(g)=n), G);

  if numelems(s) = 1 then
    return s[1]
  else
    t, r := selectremove((x -> numelems(x) = 1), s);
    if t <> [] then
      return t[1]
    else
     d := map(e -> e[2..-1] -~ e[1..-2], r):
     zip((x, y) -> if not(has(y, 1)) then x end if, r, d);
     return select(has, %, {F[]})[1]
    end if:
  end if;

end proc:

>	GeneratorsOf(10);

(1)

>	GeneratorsOf(100)

(2)

>	GeneratorsOf(142)

(3)

>	for n from 2 to 100 do printf("%4d %a\n", n, GeneratorsOf(n)) end do;

2 [2]
3 [3]

   4 [1, 3]
   5 [5]
   6 [1, 5]
   7 [2, 5]
   8 [8]
   9 [1, 8]
  10 [2, 8]
  11 [3, 8]
  12 [1, 3, 8]
  13 [5, 8]
  14 [1, 13]
  15 [2, 13]
  16 [3, 13]
  17 [1, 3, 13]
  18 [5, 13]
  19 [1, 5, 13]
  20 [2, 5, 13]
  21 [8, 13]
  22 [1, 21]
  23 [2, 21]
  24 [3, 21]
  25 [1, 3, 21]
  26 [5, 21]
  27 [1, 5, 21]
  28 [2, 5, 21]
  29 [8, 21]
  30 [1, 8, 21]
  31 [2, 8, 21]
  32 [3, 8, 21]
  33 [1, 3, 8, 21]
  34 [13, 21]
  35 [1, 34]
  36 [2, 34]
  37 [3, 34]
  38 [1, 3, 34]
  39 [5, 34]
  40 [1, 5, 34]
  41 [2, 5, 34]
  42 [8, 34]
  43 [1, 8, 34]
  44 [2, 8, 34]
  45 [3, 8, 34]
  46 [1, 3, 8, 34]
  47 [13, 34]
  48 [1, 13, 34]
  49 [2, 13, 34]
  50 [3, 13, 34]
  51 [1, 3, 13, 34]
  52 [5, 13, 34]
  53 [1, 5, 13, 34]
  54 [2, 5, 13, 34]
  55 [21, 34]
  56 [1, 55]
  57 [2, 55]
  58 [3, 55]
  59 [1, 3, 55]
  60 [5, 55]
  61 [1, 5, 55]
  62 [2, 5, 55]
  63 [8, 55]
  64 [1, 8, 55]
  65 [2, 8, 55]
  66 [3, 8, 55]
  67 [1, 3, 8, 55]
  68 [13, 55]
  69 [1, 13, 55]
  70 [2, 13, 55]
  71 [3, 13, 55]
  72 [1, 3, 13, 55]
  73 [5, 13, 55]
  74 [1, 5, 13, 55]
  75 [2, 5, 13, 55]
  76 [21, 55]
  77 [1, 21, 55]
  78 [2, 21, 55]
  79 [3, 21, 55]
  80 [1, 3, 21, 55]
  81 [5, 21, 55]
  82 [1, 5, 21, 55]
  83 [2, 5, 21, 55]
  84 [8, 21, 55]
  85 [1, 8, 21, 55]
  86 [2, 8, 21, 55]

  87 [3, 8, 21, 55]
  88 [1, 3, 8, 21, 55]
  89 [34, 55]
  90 [1, 89]
  91 [2, 89]
  92 [3, 89]
  93 [1, 3, 89]
  94 [5, 89]
  95 [1, 5, 89]
  96 [2, 5, 89]
  97 [8, 89]
  98 [1, 8, 89]
  99 [2, 8, 89]
100 [3, 8, 89]

Download fib.mw

Minor corrections...

Answered: mmcdara 7891

March 02 2024

0 0

In red in the attached file

TriagleSemblablesTest.mw

An1 := animate(display, [('f')(T2, phi)], phi = 0 .. arctan(4/3), background = T1);
...
An2 := animate(display, [('g')(T3, h)], h = 0 .. -15.8, background = T1);

Workaround...

Answered: mmcdara 7891

March 02 2024

1 0

Use intat~ instead of int :

>

restart

>	kernelopts(version)

(1)

>	with(ScientificConstants): with(Units):

>	g := evalf(Constant(g, units))

(2)

>	int(g, t*Unit('s'));

Error, (in int) integration range or variable must be specified in the second argument, got t*Units:-Unit(s)

>	# Workaround intat(g, tau=t*Unit('s')); simplify(%, 'Unit');

(3)

>	# The vector case with(VectorCalculus): SetCoordinates('cartesian[x, y, z]'):

>	a := `<,>`(0, 1, 0) *~ g

$a := Vector(3, {(1) = 0., (2) = 9.80665*Units:-Unit(m/s^2), (3) = 0.})$

(4)

>	v[0] := `<,>`(v__0, 0, 0) *~Unit('m'/'s')

v[0] := Vector(3, {(1) = `#msub(mi("v"),mi("0"))`*Unit('m'/'s'), (2) = 0, (3) = 0})

(5)

>	V := simplify(intat~(a, tau=t*Unit('s')));

dv := Vector(3, {(1) = 0., (2) = 9.80665*t*Units:-Unit(m/s), (3) = 0.})

(6)

Download Int_vector_Unit.mw

If you want to integrate twice to get the position X do this (remember that v__0 and t are dimensionless velocity and time) :

X := intat~(intat~(a, tau=t*Unit('s')), tau=t*Unit('s'))
     + 
     intat~(v[0], tau=t*Unit('s')):

simplify(X);

    v__0*t*Unit('m')*e[x] + 9.80665*t^2*Unit('m')*e[y] + 0.*e[z]

Contribution to the final answer...

Answered: mmcdara 7891

February 28 2024

3 1

I was capable to derive a surprisingly simple expression based on the hypergeometric function but not the one you want (or more precisely I got a LaguerreL based expression but the result is not correct).
In the route to the proof I have used two properties of the HermiteH polynomials (physical form) that Maple "doesn't know" and I did some transformations by hand (for instance commuting the operators Int and Sum. I don't know if this is something that can be done using directly Maple?

Maybe someone will reuse the material presented in the attached file and finalize the demonstration?
Goog luck.

>

restart;

>	with(IntegrationTools):

>	L := (m, n) -> exp(a^2)a^(abs(n-m))sqrt(2^(abs(n-m))min(n, m)!/max(n, m)!)LaguerreL(min(n, m), abs(n-m), -2*a^2);

proc (m, n) options operator, arrow; exp(a^2)*a^abs(n-m)*sqrt(2^abs(n-m)*factorial(min(n, m))/factorial(max(n, m)))*LaguerreL(min(n, m), abs(n-m), -2*a^2) end proc

(1)

>	H := 2^(-m/2-n/2)exp(lambda^2sigma^2/4)/sqrt(n!)/sqrt(m!Pi)Int(HermiteH(m,s+lambdasigma/2)HermiteH(n,s+lambdasigma/2)exp(-s^2), s=-infinity..infinity); # Set a = lambdasigma/2 H := unapply(eval(H, lambda=a2/sigma), (m, n));

2^(-(1/2)*m-(1/2)*n)*exp((1/4)*lambda^2*sigma^2)*(Int(HermiteH(m, s+(1/2)*lambda*sigma)*HermiteH(n, s+(1/2)*lambda*sigma)*exp(-s^2), s = -infinity .. infinity))/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

proc (m, n) options operator, arrow; 2^(-(1/2)*m-(1/2)*n)*exp(a^2)*(Int(HermiteH(m, s+a)*HermiteH(n, s+a)*exp(-s^2), s = -infinity .. infinity))/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2)) end proc

(2)

>	# I think this is the key property to use # # Basic Property 1: BP1m := HermiteH(m, s+a) = Sum(binomial(m, i)HermiteH(i, s)(2a)^(m-i), i=0..m); BP1n := HermiteH(n, s+a) = Sum(binomial(n, j)HermiteH(j, s)(2a)^(n-j), j=0..n);

HermiteH(m, s+a) = Sum(binomial(m, i)*HermiteH(i, s)*(2*a)^(m-i), i = 0 .. m)

HermiteH(n, s+a) = Sum(binomial(n, j)*HermiteH(j, s)*(2*a)^(n-j), j = 0 .. n)

(3)

>	# So: SBP1 := eval(H(m, n), [BP1m, BP1n]);

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*(Int((Sum(binomial(m, i)*HermiteH(i, s)*(2*a)^(m-i), i = 0 .. m))*(Sum(binomial(n, j)*HermiteH(j, s)*(2*a)^(n-j), j = 0 .. n))*exp(-s^2), s = -infinity .. infinity))/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

(4)

>

# Commute Sum operators and the Int operator

OutInt := mul(remove(has, [op(SBP1)], Int)[]);
GetI := GetIntegrand(SBP1):
SumTerms := op([1, 1], GetI)*op([2, 1], GetI):

S, R := selectremove(has, [op(%)], HermiteH):

SBP1 := OutInt * Sum(Sum(mul(R[]) * Int(mul(S[])*exp(-s^2), s=-infinity..+infinity), op([1, 2], GetI)), op([2, 2], GetI));

2^(-(1/2)*m-(1/2)*n)*exp(a^2)/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*(Sum(Sum(binomial(m, i)*(2*a)^(m-i)*binomial(n, j)*(2*a)^(n-j)*(Int(HermiteH(i, s)*HermiteH(j, s)*exp(-s^2), s = -infinity .. infinity)), i = 0 .. m), j = 0 .. n))/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

(5)

>

# If we expand SSum we will get (m+1)*(n+1) terms of the form
# T[i, j] = int(C[i, j]*HermiteH(i, x)*Hermite(j, x)*exp(-x^2), x=-infinity..+infinity).
#
# But we know that
# 1) T[i, j<>i] = 0
# 2) T[i, j=i] = sqrt(Pi)*2^i*i!
#
# It is easy to see that all terms of the form T[i, i] are got for i = 0.. min(m, n)
# Then the product of the two dpuble sums reduces to

# Basic Property 2:

BP2 := Int(HermiteH(i, s)*HermiteH(j, s)*exp(-s^2), s = -infinity .. infinity) = piecewise(j=i, sqrt(Pi)*2^i*i!, 0);

BP2 := Int(HermiteH(i, s)*HermiteH(j, s)*exp(-s^2), s = -infinity .. infinity) = piecewise(j = i, sqrt(Pi)*2^i*factorial(i), 0)

(6)

>	# Use basic property 2 into SBP1 SBP2 := eval(SBP1, BP2); # Symplify by hand SBP2 := OutInt * Sum(binomial(m, i)binomial(n, i)(2a)^(m-i)(2a)^(n-i)sqrt(Pi)2^ii!, i=0..min(m, n));

SBP2 := 2^(-(1/2)*m-(1/2)*n)*exp(a^2)*(Sum(Sum(binomial(m, i)*(2*a)^(m-i)*binomial(n, j)*(2*a)^(n-j)*piecewise(j = i, sqrt(Pi)*2^i*factorial(i), 0), i = 0 .. m), j = 0 .. n))/(sqrt(factorial(n))*sqrt(factorial(m)*Pi))

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*(Sum(binomial(m, i)*binomial(n, i)*(2*a)^(m-i)*(2*a)^(n-i)*Pi^(1/2)*2^i*factorial(i), i = 0 .. min(m, n)))/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

(7)

>	# Ask a little help (a huge one in fact) to Maple: eval(SBP2, Sum=sum);

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*Pi^(1/2)*((2*a)^(m+n)*hypergeom([-m, -n], [], (1/2)/a^2)-2*binomial(m, min(m, n)+1)*binomial(n, min(m, n)+1)*(2*a)^(m-2-2*min(m, n)+n)*factorial(min(m, n)+1)*hypergeom([1, -m+1+min(m, n), -n+1+min(m, n)], [2+min(m, n)], (1/2)/a^2)*2^min(m, n))/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

(8)

>	# Two a priori cases... which appear to give the same result mLEn := eval(SBP2, Sum=sum) assuming m <= n; mGTn := eval(SBP2, Sum=sum) assuming m > n; simplify(mLEn - mGTn);

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*Pi^(1/2)*(2*a)^(m+n)*hypergeom([-n, -m], [], (1/2)/a^2)/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

(9)

>	# Define the expression of H(m, n) this way Hermite2Hypergeom := mLEn

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*Pi^(1/2)*(2*a)^(m+n)*hypergeom([-m, -n], [], (1/2)/a^2)/(factorial(n)^(1/2)*(factorial(m)*Pi)^(1/2))

(10)

>	# Check case m < n simplify(L(1, 3), 'LaguerreL'); simplify(eval(H(1, 3), Int=int)); simplify(eval(Hermite2Hypergeom, [m=1, n=3]))

(11)

>	# Check case m = n normal(simplify(L(3, 3), 'LaguerreL')); simplify(eval(H(3, 3), Int=int)); simplify(eval(Hermite2Hypergeom, [m=3, n=3]))

(12)

>	# Check case m > n simplify(L(3, 1), 'LaguerreL'); simplify(eval(H(3, 1), Int=int)); simplify(eval(Hermite2Hypergeom, [m=3, n=1]))

(13)

>	# Now let's start from Hermite2Hypergeom and try to convert it as LaguerreL. # # Roughly: convert(Hermite2Hypergeom, `1F1`) assuming m::posint, n::posint: convert~(%, LaguerreL): Hermite2Laguerre := convert~(%, factorial);

2^(-(1/2)*m-(1/2)*n)*exp(a^2)*Pi^(1/2)*(2*a)^(m+n)*factorial(-n+m-1)*LaguerreL(m, n-m, -2*a^2)*factorial(m)*factorial(n-m)/(factorial(n)^(3/2)*(factorial(m)*Pi)^(1/2)*factorial(-n-1)*(-2*a^2)^m)

(14)

>

# That seems pretty, but there is a main problem: (-n-1)! = 0 for any positive integer n.
#
# (Note that the above expression does not account for the sign of m-n, which seems already
# quite suspect... but Hermite2Hypergeom did not neither andtheresults where correct.)

eval(Hermite2Laguerre, m=3);
eval(%, n=1);

-(1/8)*2^(-3/2-(1/2)*n)*exp(a^2)*6^(1/2)*(2*a)^(3+n)*factorial(-n+2)*LaguerreL(3, n-3, -2*a^2)*factorial(n-3)/(factorial(n)^(3/2)*factorial(-n-1)*a^6)

Error, numeric exception: division by zero

>	# I have spent a lot of time to try and understand how to get rid of this problem without any success. # # So I feel that the formula L(m, n) is not that far but I am not capable to get it. # I hope someone will and that and my work will help.

Download Q1_mmcdara.mw

Semantic correction and solution...

Answered: mmcdara 7891

February 26 2024

1 2

Writting "In the above integral, r and sigma are the random variables" is incorrect: R and Sigma are txo random variables whose realizations noted r and sigma.
Assuming that this is what you wanted to say, here is the solution

>

restart:

>	with(Statistics):

>	R := RandomVariable(Normal(k*sigma, sigma)): f__R := unapply(PDF(R, r), r);

proc (r) options operator, arrow; (1/2)*2^(1/2)*exp(-(1/2)*(-k*sigma+r)^2/sigma^2)/(Pi^(1/2)*sigma) end proc

(1)

>	Sigma := RandomVariable(LogNormal(m__Sigma, s__Sigma)): f__Sigma := unapply( (PDF(Sigma, sigma) assuming sigma > 0), sigma);

proc (sigma) options operator, arrow; (1/2)*2^(1/2)*exp(-(1/2)*(ln(sigma)-m__Sigma)^2/s__Sigma^2)/(sigma*s__Sigma*Pi^(1/2)) end proc

(2)

>	# Compute the inner integral (it could be done head-on) InnerIntegral := Int(r*f__R(r), r=-infinity..+infinity); InnerIntegral := Mean(R);

Int((1/2)*r*2^(1/2)*exp(-(1/2)*(-k*sigma+r)^2/sigma^2)/(Pi^(1/2)*sigma), r = -infinity .. infinity)

(3)

>	OuterIntegral := upperbound -> Int(sigmaf__Sigma(sigma)InnerIntegral, sigma=0..upperbound); value(OuterIntegral(+infinity)); eval(%) assuming s__Sigma > 0;

OuterIntegral := proc (upperbound) options operator, arrow; Int(sigma*`#msub(mi("f"),mi("Σ",fontstyle = "normal"))`(sigma)*InnerIntegral, sigma = 0 .. upperbound) end proc

(1/2)*2^(1/2)*k*piecewise(csgn(1/s__Sigma^2) = 1, exp(2*s__Sigma^2+2*m__Sigma)*2^(1/2)*csgn(1/s__Sigma)*s__Sigma*Pi^(1/2), infinity)/(s__Sigma*Pi^(1/2))

(4)

>	# A simpler way to get the previous result: Moment(Sigma, 2)*coeff(InnerIntegral, sigma): simplify(%); # (here the call to Statistics:-Moment automatically accounts for the conditions # on the parameters of Sigma:) print(): (attributes(Sigma)[3]):-Conditions;

(5)

>	# What happens if the integration is done for upperbound U < +infinity? value(OuterIntegral(U)): solution := simplify( (eval(%) assuming s__Sigma > 0, U > 0) )

(1/2)*k*exp(2*s__Sigma^2+2*m__Sigma)*(1+erf((1/2)*2^(1/2)*(-2*s__Sigma^2+ln(U)-m__Sigma)/s__Sigma))

(6)

>	# Now give to these parameters the values you want indets(solution, name) minus {Pi}

${U, k, m__Sigma, s__Sigma}$

(7)

>

Download integral.mw

Now, why is this semantic point important?

If r and sigma were INDEED random variables, then the inner integral (and even the outer) would make no sense. Even the simpler int(r, r=-infinity..+infinity) makes no sense.
If Sigma (I prefer using a capital do refer to a random variable) was a random variable, then
R, defined by PDF(R, r) = Normal(k*sigma, sigma) where sigma is a realization of Sigma, would not be a random variable but a random process.
A (1D) random proces is a (finite or not) collection of random variables indexed by a number (which can be itself the realization of a random variable).
Thus, "R conditionned by Sigma=sigma", which I am going to note R[sigma] is a random variable (a gaussian one with mean k*sigma and standard deviation sigma), but not R which is an infinite collection of random variables R[sigma] with sigma=0..+infinity.

Classic problems related to use 2D input...

Answered: mmcdara 7891

February 25 2024

2 2

This kind of issues are widely adressed (and answered here).
To acoid them be EXTREMELY careful when using 2D input mode in document style ... or avoid it (my advice) and use 1D input mode and worksheet instead.

The explanation is in bold brown courier font in the attached file:
(I changed the construction of U, V and Upsilon for my Maple 2015 doesn't support yours... but this is anecdotal).

>

restart

>

(1)

>

(2)

>

Phi := proc (z) options operator, arrow; z*(1-(4/3)*z^2/h^2) end proc

proc (z) options operator, arrow; z*(1-(4/3)*z^2/h^2) end proc

(3)

>

(4)

>

Q[113] := int(Q[11]*Phi(z), z = -(1/2)*h .. (1/2)*h); 1; Q[114] := int(Q[11]*Phi(z)^2, z = -(1/2)*h .. (1/2)*h); 1; Q[115] := int(Q[11]*z*Phi(z), z = -(1/2)*h .. (1/2)*h)

(5)

>

Q[553] := int(Q[55]*Phi(z), z = -(1/2)*h .. (1/2)*h); 1; Q[554] := int(Q[55]*Phi(z)^2, z = -(1/2)*h .. (1/2)*h); 1; Q[555] := int(Q[55]*z*Phi(z), z = -(1/2)*h .. (1/2)*h)

(6)

>

Q[556] := int(Q[55]*(diff(Phi(z), z)), z = -(1/2)*h .. (1/2)*h); 1; Q[557] := int(Q[55]*z*(diff(Phi(z), z)), z = -(1/2)*h .. (1/2)*h); 1; Q[558] := int(Q[55]*Phi(z)*(diff(Phi(z), z)), z = -(1/2)*h .. (1/2)*h); 1; Q[559] := int(Q[55]*(diff(Phi(z), z))^2, z = -(1/2)*h .. (1/2)*h)

(7)

>

# WHAT DID YOU DEFINE?
#
# To see really what the variables you defined are do this

lprint(anames(user));

# Observe there are 3 "Q's":
#
# on one side Q__11, Q__55
# on the other one Q
#
# This means that what appears as Q subscriped 110, 550, 113, 553, 556, ... are the same variable Q
# with indices 110, 550, 113, 553, 556, ....
#
# To see this:

indices(Q);

# Q is a table:

eval(Q);

G, h, l, Phi, rho, Q__55, W, Upsilon, U, E, m, nu, Q, L, R, mu, Q__11

table( [( 111 ) = 0., ( 550 ) = 91826.08695, ( 110 ) = 296213.1837, ( 551 ) = 0., ( 557 ) = 0., ( 556 ) = 0.1173333333e-3*Q[55], ( 559 ) = 0.9386666667e-4*Q[55], ( 558 ) = 0., ( 553 ) = 0., ( 552 ) = 0.2370337391e-3, ( 555 ) = 0.3634517333e-12*Q[55], ( 554 ) = 0.2942228318e-12*Q[55], ( 114 ) = 0.2942228318e-12*Q[11], ( 115 ) = 0.3634517333e-12*Q[11], ( 112 ) = 0.7646249649e-3, ( 113 ) = 0. ] )

(8)

>	# Note that neither [11] nor [55] indices of Q

>

beta[1] := int(mu*z^0, z = -(1/2)*h .. (1/2)*h); 1; beta[2] := int(mu*z, z = -(1/2)*h .. (1/2)*h); 1; beta[3] := int(mu*z^2, z = -(1/2)*h .. (1/2)*h)

(9)

>

beta[4] := int(mu*Phi(z), z = -(1/2)*h .. (1/2)*h); 1; beta[5] := int(mu*Phi(z)^2, z = -(1/2)*h .. (1/2)*h); 1; beta[6] := int(mu*z*Phi(z), z = -(1/2)*h .. (1/2)*h); 1; beta[7] := int(mu*(diff(Phi(z), z)), z = -(1/2)*h .. (1/2)*h)

(10)

>

beta[8] := int(mu*z*(diff(Phi(z), z)), z = -(1/2)*h .. (1/2)*h);

(11)

>

beta[11] := int(mu*(diff(Phi(z), z, z)), z = -(1/2)*h .. (1/2)*h);

(12)

>

(13)

>

I__0 := int(rho, z = -(1/2)*h .. (1/2)*h); 1; I__1 := int(rho*z, z = -(1/2)*h .. (1/2)*h); 1; I__2 := int(rho*z^2, z = -(1/2)*h .. (1/2)*h); 1; I__3 := int(rho*Phi(z), z = -(1/2)*h .. (1/2)*h); 1; I__4 := int(rho*Phi(z)^2, z = -(1/2)*h .. (1/2)*h); 1; I__5 := int(rho*z*Phi(z), z = -(1/2)*h .. (1/2)*h)

(14)

>

$eq1 := -Q__110*(diff(U(x, t), x, x))+Q__111*(diff(W(x, t), x, x, x))-Q__113*(diff(Upsilon(x, t), x, x))-Q__110*(diff(W(x, t), x))/R-`#mi("\`Q__556\`")`*Upsilon(x, t)/R+`#mrow(mi("\`Q__550\`"),mo("⋅"),mi("U"),mfenced(mrow(mi("x"),mo(","),mi("t"))))`/R^2-`#mi("\`Q__551\`")`*(diff(W(x, t), x))/R^2+`#mrow(mi("\`Q__553\`"),mo("⋅"),mi("&Upsi;",fontstyle = "normal"),mfenced(mrow(mi("x"),mo(","),mi("t"))))`/R^2+l^2*beta__1*(diff(W(x, t), x, x, x))/(4*R)-beta__7*l^2*(diff(Upsilon(x, t), x, x))/(8*R)-beta__1*l^2*(diff(U(x, t), x, x))/(8*R^2)+l^2*beta__2*(diff(W(x, t), x, x, x))/(8*R^2)-beta__4*l^2*(diff(Upsilon(x, t), x, x))/(8*R^2)+I__0*(diff(U(x, t), t, t))-I__1*(diff(W(x, t), x, t, t))+I__3*(diff(Upsilon(x, t), t, t))$

$3186210.099*Q__110*U0*exp(I*omega*t)*cos(1784.995826*x)-5687371727.*Q__111*W0*exp(I*omega*t)*cos(1784.995826*x)+3186210.099*Q__113*Upsilon0*exp(I*omega*t)*cos(1784.995826*x)-0.1784995826e-6*Q__110*W0*exp(I*omega*t)*cos(1784.995826*x)-(1/10000000000)*`#msub(mi("\`Q"),mi("556\`"))`*Upsilon0*exp(I*omega*t)*cos(1784.995826*x)+(1/100000000000000000000)*`#mrow(msub(mi("\`Q"),mi("550\`")),mo("⋅"),mi("U"),mfenced(mrow(mi("x"),mo(","),mi("t"))))`-0.1784995826e-16*`#msub(mi("\`Q"),mi("551\`"))`*W0*exp(I*omega*t)*cos(1784.995826*x)+(1/100000000000000000000)*`#mrow(msub(mi("\`Q"),mi("553\`")),mo("⋅"),mi("&Upsi;",fontstyle = "normal"),mfenced(mrow(mi("x"),mo(","),mi("t"))))`-0.4404300665e-10*beta__1*W0*exp(I*omega*t)*cos(1784.995826*x)+0.1233700550e-13*beta__7*Upsilon0*exp(I*omega*t)*cos(1784.995826*x)+0.1233700550e-23*beta__1*U0*exp(I*omega*t)*cos(1784.995826*x)-0.2202150332e-20*beta__2*W0*exp(I*omega*t)*cos(1784.995826*x)+0.1233700550e-23*beta__4*Upsilon0*exp(I*omega*t)*cos(1784.995826*x)-.2147200000*U0*omega^2*exp(I*omega*t)*cos(1784.995826*x)$

(15)

>	# WHAT ARE THE INDETERMINATES OF eq1? lprint~(indets(eq1, name));

Q__110

Q__111
Q__113
U0
W0
beta__1
beta__2
beta__4
beta__7
omega
t
x
`#mi("\`Q__551\`")`
`#mi("\`Q__556\`")`
Upsilon0
`#mrow(mi("\`Q__550\`"),mo("⋅"),mi("U"),mfenced(mrow(mi("x"),mo(","),mi("t"))))`
`#mrow(mi("\`Q__553\`"),mo("⋅"),mi("&Upsi;",fontstyle = "normal"),mfenced(mrow(mi("x"),mo(","),mi("t"))))`

(16)

>	# Compare to the table Q whose indices are remembered here {indices(Q, nolist)}

(17)

>	# So the answer to your question is obvious: # # You defined Q[110] and use Q__110 in eq1 (same thing for the other quantities) # # And the remedy is simple: be more carefull when you se 2D mode # And a good advice (IMO), let's drop it and use 1D input mode instead to avoid this kind of issues.

Download TRASH_mmcdara.mw

Note that your problem does not limit tho the Qs but extends to the betas.

The problem with the 2D inputs is that what you see is not what you get.
Here you falsely believed that the variables in eq1 where those you previously defined.

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

7891 Reputation

22 Badges

MaplePrimes Activity

These are answers submitted by mmcdara

Here is a partial result...

A thirteen steps solution......

A solution...

Explanation...

If I may offer some advice...

gamma l!!!...

Why did you delete your previous questio...

A way, likely among many...

Confusion between variance and raw momen...

If I am not mistaken...

Minor corrections...

Workaround...

Contribution to the final answer...

Semantic correction and solution...

Classic problems related to use 2D input...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again