Integral of a gaussian function wrong answer

I am trying to calculate the integral of the following Gaussian function:

syms x A
assume(x,'real');
assume(A,'real');
assumeAlso(A>0);
% Expression to integrate
eq = x^14*exp(-x^2/2)*exp(-A*x);
% Integration
ans = int(eq,x,-inf,inf)

Matlab gives me an answer containing hypergeo and erf. However, there is a regular algebraic solution without the needs of hypergeo or erf. (For example type "integral of (x^14*exp(-x^2/2)*exp(-A*x)) from -inf to inf" into wolfram alpha to get algebraic solution)

Is there a way to get to the algebraic answer in Matlab?

Answer 1

[Karan Gill]

I could get the algebraic answer but it required calling "simplify" twice. Need to investigate this further.

>> syms x A
assume(x,'real');
assume(A,'real');
assumeAlso(A>0);
% Expression to integrate
eq = x^14*exp(-x^2/2)*exp(-A*x);
sol = int(eq,x,-Inf,Inf);    % long expression

>> sol1 = simplify(sol,'Steps',500); % first simplify step

>> sol2 = simplify(sol1,'Steps',200) % second simplify step
sol2 =
2^(1/2)*pi^(1/2)*exp(A^2/2)*(A^14 + 91*A^12 + 3003*A^10 + 45045*A^8 + 315315*A^6 + 945945*A^4 + 945945*A^2 + 135135)

[Matthias Freiherr von Andrian-Werburg]

Thank you! I did not experiment with the "step" option of simplify before.

[Karan Gill]

You can find the option in 4th example under simplify here. If the example is hard to find or read, that would be useful to know. Unfortunately, sometimes these options do get lost in longer doc pages.

Thanks,

Karan (Symbolic doc)

[Matthias Freiherr von Andrian-Werburg]

The "step" option is clearly documented, as you point out. I just did not think of using it, because so far I was able to simplify the expressions I was dealing with with the default "step" option.

Answer 2

[Matthias Freiherr von Andrian-Werburg]

Follow up:

The provided solution does not work, when I multiply the expression with an exponent of a constant:

syms x A1 A2 B C
assume(x,'real');
assume(A1,'real');
assume(A2,'real');
assume(B,'real');
assumeAlso(A1>0);
assumeAlso(A2>0);

% Expression to integrate
eq0 = x^14*exp(-A2*x)*exp(-x^2/2); % Simplest expression
eq1 = B*eq0; % Multiplying with constant
eq2 = exp(-A1)*eq0; % Multiplying with exp(constant)
eq3 = B*exp(-A1)*eq0; % Multiplying with both

% Integration (comment out all lines, but one)
% sol1 = int(eq0,x,-inf,inf); % -> works
% sol1 = int(eq1,x,-inf,inf); % -> works
% sol1 = int(eq2,x,-inf,inf); % -> does not work
sol1 = int(eq3,x,-inf,inf); % -> does not work

% Simplification to get from hypergeom and erf solution to algebraic solution
sol2 = simplify(sol1,'steps',1000)
sol3 = simplify(sol2,'steps',200)

Is there a quick simplification to get to the algebraic solution?

Alternatively, I could imagine splitting up the integrand into independent and dependent parts of x and only integrating the dependent parts. Is there a way to use the "children" function to split expressions into dependent and independent parts?

[Matthias Freiherr von Andrian-Werburg]

In case anybody ever runs into the same problem, here is how I worked around the issue:

syms x A1 A2 B C
assume(x,'real');
assume(A1,'real');
assume(A2,'real');
assume(B,'real');
assumeAlso(A1>0);
assumeAlso(A2>0);

% Expression to integrate
eq0 = x^14*exp(-A2*x)*exp(-x^2/2); % Simplest expression
eq1 = B*exp(-A1)*eq0; % Multiplying with constant and exponent of constant

% Integration
sol1 = int(eq1,x,-inf,inf); % -> solution with hypergeom and erf

% Simplification to get from hypergeom and erf solution to algebraic solution
sol2 = simplify(sol1,'steps',1000)
sol3 = simplify(sol2,'steps',200)
% -> does not work

% Factor expression
fact = factor(eq1);
% find factors with x and without x
fact_x = 1;
fact_no_x = 1;
for ii = 1:length(fact)
    if ~isempty(find(symvar(fact(ii))==x))
        fact_x = fact_x * fact(ii);
    else
        fact_no_x = fact_no_x * fact(ii);
    end
end

% Integrate factorized expression with x
sol4 = int(fact_x,x,-inf,inf); % -> solution with hypergeom and erf

% Simplification to get from hypergeom and erf solution to algebraic solution
sol5 = simplify(sol4,'steps',1000)
sol6 = simplify(sol5 * fact_no_x,'steps',200) % Bring back factors without x
% -> works