Explicit integral could not be found.

Hello,

I am trying to get a solution for below problem, which is to evaluate double integrals of function f(x,y) where the first integral boundary is in variable y. Though the explicit integral could not be found, I expect numerical answer still exists; however I don't know how to get that. I tried double and VPA functions but they didn't work.

___________________________________________________________________________

syms x y

f = x^3*exp(1/x^2)*(x^(1/2)/2)^(1/2);

F = int(f,x,0,y);

ans = int(F,y,0,100); ___________________________________________________________________________

Warning: Explicit integral could not be found.

Warning: Explicit integral could not be found.

ans = int(int(x^3*exp(1/x^2)*(x^(1/2)/2)^(1/2), x = 0..y), y = 0..100)

___________________________________________________________________________

 VPA(ans)
ans =
numeric::int(numeric::int(x^3*exp(1/x^2)*(x^(1/2)/2)^(1/2), x = 0..y), y = 0..100)

Answer 1

[Babak]

Mathematica finds the indefinite integral of f over x

http://integrals.wolfram.com/index.jsp?expr=x%5E3*exp%281%2Fx%5E2%29*%28x%5E%281%2F2%29%2F2%29%5E%281%2F2%29&random=false

Take this and evaluate it for the boundaries (i.e. x=y and x=0) and subtract the results. this will give you the definite integral from x=0 to x=y then again try to integrate the result over y.

[Babak]

To evaluate the value of int(f) at x=0 you need to compute limit(int(f)) when x->0+, because I can see that int(f) is singular at x=0.

[David Yoo]

Thank you for the answer. However, is there any way to solve this kind of problem in Matlab?

[Babak]

Try Matt's solution.

Answer 2

[Matt Fig]

Are you sure about those limits?

Here is how one would ordinarily solve such a problem numerically. Note I use limits 1-y and 2-10:

fh = @(y) quadl(@(x) x.^3.*exp(1./x.^2).*(x.^(1/2)/2).^(1/2),1,y);
quadv(fh,2,10)

.

.

.

Let's just check the method with an easy problem. Say we want to solve:

int(int(x,0,y),1/2,1)

By simple calculus this is:

int(y^2/2,1/2,1) = 7/48

So in MATLAB

fh = @(y) quadl(@(x) x,0,y);
quadv(fh,1/2,1) % Equal to 7/48 as expected...