matlabnumerical-integrationnumerical-computing

How to solve a double integral with any code?

I have an double integral:

$$\int_0^{\infty} y^2 e^{\int_0^y e^{-x^2} d x} d y$$

Because exp(−𝑥^2) is a non-integrable function, I have tried to solve this using the quadgk function in MATLAB, but I don't get a good result. Changing the integral's upper limit from infinity to some exact value may be a good compromise.

I had fitted

$${\int_0^y e^{-x^2} d x}$$

with a polynomial fitting, so the whole formula can be analytic. However, sometimes the polynomial fitting is badly conditioned, so I need a better idea to get a better solution.

Solution

  • You can solve this numerically using two calls to integral like this:

    g=@(x)exp(-x.^2);
h=@(y)y.^2;
a=20;
z=integral(@(y)h(y).*exp(integral(g,0,y)),0,a,'ArrayValued',true);

    You need to use the 'ArrayValued' option on the outer call to integral to force the function to run in a non-vectorized mode so the underlying algorithm doesn't try to pass a vector argument in which would be invalid for the upper integration limit of the inner integral.

    And you can verify this result symbolically using Matlab's Symbolic Math toolbox and either vpa or double to evaluate the integral:

    syms x y real;
g=exp(-x^2);
h=y^2;
a=20;

z=int(h*exp(int(g,x,0,y)),y,0,a)
vpa(z)    % evaluate and output as variable precision arithmetic
double(z) % evaluate and output as floating point