How to solve a double integral with any code?

Question

I have an double integral:

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

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.

Answer 1

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

Answer 2

without using any int function,this solve the double integral from zero to a(here a equals 20);

clear;clc;
N = 1e7;
a=20;
h = a/N;
x = 0:h:a;
y = x;

z = zeros(size(x));
for i = 2:N+1
    z(i) = z(i-1) + exp(-y(i-1)^2)*h;
end
%parameter z(end) means int(exp(-x^2),0,a) 

%the integrand means the bulk of integrand which can be express  
%any uplimit below a;
integrand = y.^2 .* exp(z);

%It just like the definite integral meanning expresses area as a 
%combination of infinitely many vertically-oriented rectangles. 
result= h *sum(integrand(1:end));
disp(result);