Numerical Integration Gives 0

Hello everyone,

I'm implementing the model proposed by Apsley (1975) Paper Tile: "Temperature- and field-dependence of hopping conduction in disordered systems, II"

In this I need to numerically integrate some quantities, which happen to need a double numerical integration involving the fermi probability function and interdependent integral limits

I try to define them in the following code (fun 1 through 4, I1 through 4)

I understand the mathematics are not as straightforward but I would appreciate some feedback on if my code has some implementation errors!

I get I1 through 4 = 0 for multiple cases which is unphysical and warnings on I4.

Any feedback is welcome!

Thank you!

---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

for j = 0:1:N

    F = Fmin + j*Fstep; % Update the field values

    bet = (F*e) / (2*a*Uth);
    %   Use Arsley (13), (14)

    P = (1 + bet/2) / (1 + bet)^2;
    Q = 3*bet/2 + 1;

    for j2 = 1:1:Nstate

        Upr_loc = Upr(1,j2);

        if Upr_loc >=  0

            R4dnn = ( 2/(K*(P+Q)) )^(1/4);
            % R4dnn, Upr_loc, bet,
            fun1 = @(Vpr, theta) ( 1 - 1 ./ (1 + exp(Vpr-Efpr)) ) .* ( (R4dnn - Vpr + Upr_loc) ./ (1 + bet .* cos(theta) ) ).^3 .* sin(theta) .* cos(theta);
            fun2 = @(Vpr, theta) ( 1 - 1 ./ (1 + exp(Vpr-Efpr)) ) .* R4dnn^3 .* sin(theta) .* cos(theta);
            fun3 = @(Vpr, theta) ( 1 - 1 ./ (1 + exp(Vpr-Efpr)) ) .* ( (R4dnn-Vpr+Upr_loc) ./ (1+bet .* cos(theta) ) ).^3 .* sin(theta);
            fun4 = @(Vpr, theta) R4dnn^.2 .* sin(theta);

            % Integrals assume Constant Density of States N(Vpr) = Nt !
            % hence effectively I1....I4 represent I1/Nt. For a non
            % constant density of states it has to be taken into account!

            Vprtheta = @(theta) Upr_loc - R4dnn*bet*cos(theta); 
            Vprmax = Upr_loc + R4dnn;
            Vprinf  = -Inf;

            I1 = integral2(fun1,0,pi,Vprtheta,Vprmax); 
            I2 = integral2(fun2,0,pi,Vprinf,Vprtheta);
            I3 = integral2(fun3,0,pi,Vprtheta,Vprmax);
            I4 = integral2(fun4,0,pi,Vprinf,Vprtheta); % Eval not good. Why ?

            RFpr = (I1 + I2)/(I3 + I4);

            flag = 1;
        else
            R4dnn = ( 2/(K*(P+Q)) )^(1/4) -((P+1)*Upr)/(P+Q); 
            flag  = 0;  
        end;

    end;

    Fvec(j+1,1) = F;

end;

Answer 1

[John D'Errico]

This is difficult to answer, without understanding the mathematics behind your code.

Perhaps the VERY most common reason why a numerical integration returns a zero result when one would expect something else is a simple one. Numerical integrations are simple tools. They sample the function. If the function is zero, or essentially so at every point they sample, they give up. Hey! Its zero. I know what the integral of zero is. ZERO.

The problem is that too often somebody throws what is essentially a delta function spike into the mix, and expect the integration to survive. Zero everywhere, except for some TINY region, where it is non-zero. Yeah, what do you expect to see happen? How does the integration tool know that this problem is any different from the one I mentioned before, where the function WAS zero everywhere?

So after due diligence, an integration tool will concede defeat when it sees only zeros everywhere. It returns zero.

A solution to the delta function thing is to use waypoints in the solver, if it admits them. That is, force the tool to look at that area as a point of importance.

As I said, this is the most common "failure" mode when people work with numerical integration tools, and most especially when the word probability appears in the conversation. Do I know this to be the case? No, but it is the very first thing I would check.

If not that, there are other issues to consider. Have you written the code properly? Hey, I don't know, at least without carefully reading it.

[Tomer Hamam]

I followed this post on a similar problem and you were spot on. I opted for the usage of 'Waypoints' in the 'integral' function options to direct him to intervals which solved the issue for my problem.

[John D'Errico]

'waypoints' are a great way to help the integral tool to know where to look. As long as you know where something is happening, that will fix it.

Answer 2

[Marios Barlas]

Thanks for the answers guys!

In any event, the sin is coupled with other functions and it should not give zero in all cases even in the [0,pi] range. This theory is the one who have Dr. Mott his Nobel so I'm willing to trust that they did something right and it's my own inabiliy to properly solve it :P

I upload the integrals in symbolic form just to make it easier to understand the integrals. forget about the N(V') function in my approximation it's constant for certain reasons. the rest of the implicated parameters are constants for the integration. Or at least should be as I defined my code!

Thanks a lot anyways!