Two variables integration matlab

Question

I'm trying to solve physical problem in quantum mechanics of helium atoms, the solution require numerical integration over 2 variables. However when i'm trying to run the next code

e = 1.6e-19;
eps0 = 8.85e-12;
a0 = 0.53e-10;
c = e^2/(4*pi*eps0);
psi1 =@(r) ((2^1.5)/(sqrt(pi)*a0^1.5))*exp(-2*r/a0);
psi2 =@(r) ((2^1.5)/(4*sqrt(2*pi)*a0^1.5))*(2-r/a0)*exp(-r/a0);
I =@(r1,r2) c*(psi2(r1).*psi1(r2).*(1./abs(r1-r2)).*psi2(r1).*psi1(r2));
J =@(r1,r2) c*(psi2(r1).*psi1(r2).*(1./abs(r1-r2)).*psi1(r1).*psi2et(r2));
integral2(I,0,inf,0,inf)
integral2(J,0,inf,0,inf)

matlab saying

Error using  * 
Inner matrix dimensions must agree.

Error in qmex7>@(r)((2^1.5)/(4*sqrt(2*pi)*a0^1.5))*(2-r/a0)*exp(-r/a0)

Error in qmex7>@(r1,r2)c*(psi2(r1).*psi1(r2).*(1./abs(r1-r2)).*psi2(r1).*psi1(r2))

Error in integral2Calc>@(y)fun(xi*ones(size(y)),y) (line 18)
    @(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions), ...

Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 83)
        [q,errbnd] = vadapt(@AToInfInvTransform,interval);

Error in integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions)

Error in integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x)) (line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 83)
        [q,errbnd] = vadapt(@AToInfInvTransform,interval);

Error in integral2Calc>integral2i (line 20)
[q,errbnd] = integralCalc(innerintegral,xmin,xmax,opstruct.integralOptions);

Error in integral2Calc (line 7)
    [q,errbnd] = integral2i(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 106)
    Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

any idea why?

Answer 1

The error message is because you forgot the . before * before exp in psi2.

Once you fix that, MATLAB will be none too happy about the singularity you get dividing by abs(r1-r2) in I and J. If (1./abs(r1-r2)) is changed, for instance, to (1./(abs(r1-r2) + 1e-20)) , MATLAB will produce an answer of 0 for integral2(I,0,inf,0,inf). I'll let you figure out how to fix that properly.

You have horrible scaling to the max, by having a0 = 0.53e-10 in the denominator of the argument of exp. Hence, for instance, psi1(1e-8) = 5.39e-149, and only gets smaller as its argument increases. And this despite a multiplicative factor (2^1.5)/(sqrt(pi)*a0^1.5) = 4.14e15, which, though, gets overwhelmed by the exponential of a large negative number. All the real action takes place at even smaller argument values than 1e-8. For instance, psi1(1e-10) = 9.60e+13, psi1(1e-9) = 0.17. So I think you have some reformulation needed to get this to the point where the integral can be evaluated in a numerically sound way.