0

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.

horchler
  • 18,384
  • 4
  • 37
  • 73
Leo_ma
  • 3
  • 2
  • Please update your question to be more explicit with what you are looking for – as is, it's a bit open ended. It's also hard for users to help you without a better understanding of h(y) and g(x). Can you provide a minimal example? – horchler Mar 01 '23 at 18:31
  • Also, have you tried [`integral`](https://www.mathworks.com/help/matlab/ref/integral.html) instead of `quadgk`? – horchler Mar 01 '23 at 18:32
  • to Leo_ma, hi, why don't you try to 1st solve int(,y,g,x) and then put that solution, using parameters, into the outer integral? Also, you mention h and g are 'complicated'. How complicated can they be? Mission Impossible? complicated, impossible, problemo, danger, may-do-boom .. whatever : if you divide complicated problems into smaller simpler ones, and put enough know-how and effort, it may even all turn out to be a walk in the park. – John BG Mar 01 '23 at 23:33
  • to horchler, thanks for your comment, and I have already provide a minimal example. – Leo_ma Mar 02 '23 at 04:30
  • to john BG, now I have the feeling to be a walk in the park after divided the complicated problems. thanks for your advise. – Leo_ma Mar 02 '23 at 04:32
  • Using the Symbolic toolbox you can show that `syms x y real;` `int(exp(-x^2),x,0,y)` returns `(pi^(1/2)*erf(y))/2`. – horchler Mar 02 '23 at 14:29
  • However, you'll see that `int(y^2*int(exp(-x^2),x,0,y),y,0,Inf)` returns `Inf`. So, as specified your integral is unbounded. If you replace the upper bound of the outer integral with a free parameter, e.g., `a`, you can obtain and expression in terms of that variable: `(exp(-a^2)*(a^2 + 1))/6 + (a^3*pi^(1/2)*erf(a))/6 - 1/6`. – horchler Mar 02 '23 at 14:33
  • to john BG,I confirmed your method to solve the int(y^2*int(exp(-x^2),x,0,y),y,0,a), that's right. But you lack 'exp' before 2nd int, after adding, the integral should be int(y^2.*exp(int(exp(-x^2),x,0,y)),y,0,a),and returns int(y^2*exp((pi^(1/2)*erf(y))/2), y, 0, a). if a equals 20, the return is int(y^2*exp((pi^(1/2)*erf(y))/2), y, 0, 20). that's not an exact number. what should we do? – Leo_ma Mar 02 '23 at 16:00
  • Here is the another method to solve the double integral to get an exact number(a=20):(https://github.com/saicma/fundamental_example/blob/main/double_integral.m), return 2.593240654688026e+04.@horchler@John BG – Leo_ma Mar 02 '23 at 16:11
  • @Leo_ma sorry, missed that extra `exp`. See my answer below for how to evaluate this numerically or symbolically. – horchler Mar 02 '23 at 22:03
  • Your link returns "404 Not Found" for me - seems to be a private repo maybe? My answer below also returns a different value that you obtained for `a=20`. I get `6.4689e+03`. – horchler Mar 02 '23 at 22:10
  • sorry, my link is a private repo, you can't to access. Otherwise my old code's answer was wrong compared with your answer. so I corrected it . My answer below returns the same value that you obtained for a=20. I also get 6.4689e+03. Last but not least, your answer is so pretty and correct, it help me a lot. thanks! @horchler – Leo_ma Mar 03 '23 at 05:55

2 Answers2

0

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
horchler
  • 18,384
  • 4
  • 37
  • 73
0

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);
Leo_ma
  • 3
  • 2
  • Thank you for your answer! Please try to enrich it with some more comments on what the your code is doing. – Georg W. Mar 06 '23 at 08:23
  • Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers [in the help center](/help/how-to-answer). – Georg W. Mar 06 '23 at 08:23