2
$\begingroup$

I have the following code:

substitute = {x_ Conjugate[x_] -> Abs[x]^2};

S = Sqrt[2]/2*{{1 + Conjugate[δ], 0}, {0, 1 - Conjugate[δ]}}; 

k = (1/Sqrt[2])*{{S[[1, 1]] + S[[2, 2]]}, {S[[1, 1]] - S[[2, 2]]}, {2 S[[1, 2]]}} //
     Simplify;

T0 = Dot[k, ConjugateTranspose[k]];

R[ψ_] := {{1, 0, 0}, {0, Cos[2 ψ], Sin[2 ψ]}, {0, -Sin[2 ψ], Cos[2 ψ]}};

T[ψ_] := Dot[R[ψ], T0, Refine[ConjugateTranspose[R[ψ]], ψ ∈ Reals]] /. substitute;

p[x_, mu_, k_] := 2*PDF[VonMisesDistribution[2 mu, k], 2 x] // Simplify;

TvolNRS = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};

TvolNRS[[1, 1]] = Integrate[p[ψ, μ, κ]*T[ψ][[1, 1]], {ψ, -Pi/2, Pi/2},
  Assumptions -> -Pi/2 <= μ <= Pi/2] // FullSimplify

enter image description here

How can I force mathematica to calculate the highlighted integral and print the result explicitly as (a rather simple) function of the parameters $\mu$ and $\kappa$?

$\endgroup$
10
  • $\begingroup$ You seem to have got the result. What do you mean by "print"? $\endgroup$ Oct 10 '15 at 17:49
  • $\begingroup$ @AlexeiBoulbitch I mean to calculate the result of the integral and print it explicitly $\endgroup$ Oct 10 '15 at 20:27
  • $\begingroup$ @MichaelE2 I have edited my question and added the codes that are necessary for running $\endgroup$ Oct 10 '15 at 20:29
  • $\begingroup$ Your integral after removing constant parameters comes down to Integrate[E^Cos[x], x], which returns unevaluated, indicating that Mathematica does not know how to find the antiderivative. Do you know if there is one in terms of standard functions? $\endgroup$
    – Michael E2
    Oct 10 '15 at 20:37
  • 1
    $\begingroup$ Yes, that's what I mean. Unless somebody knows a symbolic solution.... $\endgroup$
    – Michael E2
    Oct 10 '15 at 21:02
3
$\begingroup$

Doing symbolic manipulations with Integrate isn't really possible so I'll just talk this out.

integrand = p[ψ, μ, κ]*T[ψ][[1, 1]];
Reduce[integrand[[1, 1, 2]], ψ][[2]]
(* 1/2 (-π + 2 μ) <= ψ <= 1/2 (π + 2 μ) *)

This suggests that we shift the integration variable by μ:

integrand = MapAt[Reduce[#, y] &, integrand /. ψ -> μ + y // Simplify, {1, 1, 2}]

enter image description here

This shifts the limits of integration to

# - μ & /@ {ψ, -π/2, π/2} /. ψ -> μ + y
(* {y, -(π/2) - μ, π/2 - μ}  *)

However, due to the constraints on y from the definition of the integrand, the actual limits are given by

{y, Max[-(π/2), -(π/2) - μ], Min[π/2, π/2 - μ]}

Visualizing the integration region for different choices of μ:

enter image description here


So, the problem has been mapped onto the following:

g[μ_ /; -π/2 <= μ <= π/2] := Integrate[integrand, {y, Max[-(π/2), -(π/2) - μ], Min[π/2, π/2 - μ]}]

At μ == 0, of course, the integral evaluates to 1, since

Integrate[Exp[κ Cos[2 y]], {y, -π/2, π/2}]
(* π BesselI[0, κ] *)

I know of no closed-form analytic solution for the integral. Interestingly enough, however, certain values work:

g[π/4] // Expand
(* 3/4 + StruveL[0, κ]/(4 BesselI[0, κ]) *)

The StruveL function shows up, which I've actually used before.


This all suggests that you have to do things numerically. The integral has been made simpler, so perhaps this was a help.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.