Subject: problem with integration of gaussian distribution
From: Stavros Macrakis
Date: Mon, 6 Jul 2009 15:45:23 -0400
On Mon, Jul 6, 2009 at 3:20 PM, Riemann, Robert <
robert.riemann at physik.hu-berlin.de> wrote:
> Maxima computes the integration.
> The result may be right, but is to complicated to interpret.
> Exponentialze doesnt seem to work for this solution.
>
> I thought i could help maxima with some "wise" ;) substitutions.
> I found out that Maxima solves something like this as expected:
> integrate(exp(-x^2-2*x),x,minf,inf);
>
Well, this has no parameters, while your integral has 6 (A, m, d, P_0, h,
t). So it is not surprising that the result is more complicated.
> What must I do, to let Maxima sort some summands in powers of p?
>
I'm sorry, I don't understand.
> Is it possible to force Maxima to use a given substitution to make the
> terms easier?
>
I'm sorry, I don't understand. You can do any substitutions you want -- see
? changevar. But you can't control how the 'integrate' command performs the
integral internally.
The simplest form I could find with basic Maxima functions, by the way, used
scanmap(gfactor,...). See below.
-s
-------------------------
(%i1) phi(p)*%e^(%i*(p*x-p^2*t/(2*m))/h)/(2*%pi*h);
(%o1) %e^(%i*(p*x-p^2*t/(2*m))/h-d^2*(p-p_0)^2/h^2)*A/(2*%pi*h)
(%i2) integrate(%,p,minf,inf);
Is %e^(d^2*p_0/h^2)-1 positive, negative, or zero?
p;
Is t positive, negative, or zero?
p;
Is x positive, negative, or zero?
p;
Is p_0 positive, negative, or zero?
p;
Is
sqrt(sqrt(h^2*t^2+4*d^4*m^2)-2*d^2*m)*(sqrt(h^2*t^2+4*d^4*m^2)+2*d^2*m)^(3/2)-h^2*t^2
positive, negative, or zero?
p;
(%o2) (%i*(sqrt(%pi)*h*sqrt(m)*sqrt(sqrt(h^2*t^2+4*d^4*m^2)+2*d^2*m)
*%e^(-2*d^2*m^2*x^2/(2*h^2*t^2+8*d^4*m^2)
+4*d^2*m*p_0*t*x/(2*h^2*t^2+8*d^4*m^2)
+8*d^6*m^2*p_0^2/(2*h^4*t^2+8*d^4*h^2*m^2)-d^2*p_0^2/h^2)
*sin(m*(h^2*t*x^2+8*d^4*m*p_0*x-4*d^4*p_0^2*t)
/(2*h*(h^2*t^2+4*d^4*m^2)))
/sqrt(h^2*t^2+4*d^4*m^2)
-sqrt(%pi)*h*sqrt(m)*sqrt(sqrt(h^2*t^2+4*d^4*m^2)-2*d^2*m)
*%e^(-2*d^2*m^2*x^2/(2*h^2*t^2+8*d^4*m^2)
+4*d^2*m*p_0*t*x/(2*h^2*t^2+8*d^4*m^2)
+8*d^6*m^2*p_0^2/(2*h^4*t^2+8*d^4*h^2*m^2)-d^2*p_0^2/h^2)
*cos(m*(h^2*t*x^2+8*d^4*m*p_0*x-4*d^4*p_0^2*t)
/(2*h*(h^2*t^2+4*d^4*m^2)))
/sqrt(h^2*t^2+4*d^4*m^2))
+sqrt(%pi)*h*sqrt(m)*sqrt(sqrt(h^2*t^2+4*d^4*m^2)-2*d^2*m)
*%e^(-2*d^2*m^2*x^2/(2*h^2*t^2+8*d^4*m^2)
+4*d^2*m*p_0*t*x/(2*h^2*t^2+8*d^4*m^2)
+8*d^6*m^2*p_0^2/(2*h^4*t^2+8*d^4*h^2*m^2)-d^2*p_0^2/h^2)
*sin(m*(h^2*t*x^2+8*d^4*m*p_0*x-4*d^4*p_0^2*t)/(2*h*(h^2*t^2+4*d^4*m^2)))
/sqrt(h^2*t^2+4*d^4*m^2)
+sqrt(%pi)*h*sqrt(m)*sqrt(sqrt(h^2*t^2+4*d^4*m^2)+2*d^2*m)
*%e^(-2*d^2*m^2*x^2/(2*h^2*t^2+8*d^4*m^2)
+4*d^2*m*p_0*t*x/(2*h^2*t^2+8*d^4*m^2)
+8*d^6*m^2*p_0^2/(2*h^4*t^2+8*d^4*h^2*m^2)-d^2*p_0^2/h^2)
*cos(m*(h^2*t*x^2+8*d^4*m*p_0*x-4*d^4*p_0^2*t)/(2*h*(h^2*t^2+4*d^4*m^2)))
/sqrt(h^2*t^2+4*d^4*m^2))
*A
/(2*%pi*h)
(%i3) scanmap(gfactor,exponentialize(%));
(%o3) sqrt(m)*sqrt((h*t-2*%i*d^2*m)*(h*t+2*%i*d^2*m))
*(sqrt(sqrt((h*t-2*%i*d^2*m)*(h*t+2*%i*d^2*m))+2*d^2*m)
-%i*sqrt(sqrt((h*t-2*%i*d^2*m)*(h*t+2*%i*d^2*m))-2*d^2*m))
*%e^(%i*(h*m*x^2-4*%i*d^2*m*p_0*x+2*%i*d^2*p_0^2*t)/(2*h*(h*t-2*%i*d^2*m)))*A
/(2*sqrt(%pi)*(h*t-2*%i*d^2*m)*(h*t+2*%i*d^2*m))
(%i4)