2
$\begingroup$

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?

$\endgroup$
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.