compute integral with Maxima



thanks to dan for the reference.

the problem is then
g(x) := exp( -( (a-b*cos(x))^2 + (c-b*sin(x))^2 ) );
integrate((g(x))^2, x, 0, 2*%pi);
or
integrate(exp( - 2* (a-b*cos(x))^2 -2*(c-b*sin(x))^2 ),x,0,2*%pi );

with the answer
2*%pi*%e^(-2*c^2-2*b^2-2*a^2)*I0(4*b*sqrt(c^2+a^2)),

with I0(x) the same as bessel_i(0,x).

A solution which avoids looking up trig identities, but
has to use exp(z*cos(th)) expansion in terms of
I(k,z) from Abramowitz and Stegun 9.6.34, page 376,
(Dan's reference):

(%i1) display2d:false$

arge is argument of the exponential in the integrand

(%i2) arge:(-2*(c-b*sin(x))^2 -2*(a-b*cos(x))^2 ,trigsimp(%%));
(%o2) 4*b*c*sin(x)+4*a*b*cos(x)-2*c^2-2*b^2-2*a^2

exp(arge1) comes outside the integral
(%i3) arge1:rest(arge,2);
(%o3) -2*c^2-2*b^2-2*a^2

the other factor is exp(arge2) which remains to
be integrated

(%i4) arge2:rest(arge,-3);
(%o4) 4*b*c*sin(x)+4*a*b*cos(x)
the trig identity replaces arge2 by u*cos(x-v)

(%i5) eqns:arge2 = u*cos(x-v);
(%o5) 4*b*c*sin(x)+4*a*b*cos(x) = u*cos(x-v)

since cos(x) and sin(x) coefficients must be zero
this is actually two equations

(%i6) eq1:integrate(eqns*sin(x),x,0,2*%pi)/(u*%pi);
(%o6) 4*b*c/u = sin(v)

(%i7) eq2:integrate(eqns*cos(x),x,0,2*%pi)/(u*%pi);
(%o7) 4*a*b/u = cos(v)

solve barfs
(%i8) solve([eq1,eq2],[u,v]);
(%o8) []

so do it "by hand"

(%i9) eq1/eq2;
(%o9) c/a = sin(v)/cos(v)
(%i10) vsoln: ratsubst(tan(v),sin(v)/cos(v),%);
(%o10) c/a = tan(v)

(%i11) trigsimp(eq1^2 + eq2^2);
(%o11) (16*b^2*c^2+16*a^2*b^2)/u^2 = 1
(%i12) solve(%,u);
(%o12) [u = -4*b*sqrt(c^2+a^2),u = 4*b*sqrt(c^2+a^2)]
(%i13) usoln:second(%);
(%o13) u = 4*b*sqrt(c^2+a^2)

noun form solution for the factor we still need to integrate:

(%i14) i1:'integrate( I0(u)+2*'sum(I(k,u)*cos(k*x-k*v),k,1,inf),x,0,2*%pi)$
(%i15) (declare(k,integer),assume(k>0))$

(%i16) i1,integrate;
(%o16) 2*%pi*I0(u)
(%i17) %,usoln;
(%o17) 2*%pi*I0(4*b*sqrt(c^2+a^2))

mult. by factor pulled outside to get the
whole integral
(%i18) isoln:exp(arge1)*%;
(%o18) 2*%pi*%e^(-2*c^2-2*b^2-2*a^2)*I0(4*b*sqrt(c^2+a^2))

check numerically at one set of (a,b,c)

(%i19) val1:[a=1/sqrt(2),b=1/2,c=1/sqrt(2)]$
(%i20) subst(val1,isoln);
(%o20) 2*%e^-(5/2)*%pi*I0(2)
(%i21) ratprint:false$
(%i22) ratsubst(bessel_i(0,2.0),I0(2),%o20);
(%o22) 12752*%e^-(5/2)*%pi/2797
(%i23) float(%);
(%o23) 1.175708087416373

now do the integral numerically from the
start

(%i24) subst(val1,arge);
(%o24) 2*sin(x)/sqrt(2)+4*2^-(3/2)*cos(x)-5/2
(%i25) farge:float(%);
(%o25) 1.414213562373095*sin(x)+1.414213562373095*cos(x)-2.5
(%i26) quad_qags(exp(farge),x,0,2*%pi);
(%o26) [1.175708104128904,7.3427092520870032E-13,105,0]
so numerical agreement.
 
Ted Woollett