quad2d problem

Hi Folks I am getting used to the quad2d function on matlab. I want to pefrom a double integral and I choose a non trivial case which is the ionization of hydrogen 1s electron by a fast proton based on the first Born approximation. The integral is well known and can be find at McDowell and Colleman book "introduction to the theory of ion-atom collision", Eq 7.2.32. My code is very short and works for energies up to 50keV. Above it the Integral reaches the maximum number of function evaluations (and the result fails the global error test) and the values are not correct. The only thing that changes from one energy to another is the lower limit integrand which dpends on de velocity consequently with the energy. I tried to comment as much as possible Thank you All units are in atomic unit = cross section are in (Pi*ao^2)

   function integral_dupla_v2
  fclose all; close all; clear all;
  %Ent = [20 25 35 50 70 80 100 120 140 160 180 200 225 300 500 600 800 1000 2000 4000];
    Ent = [1 2 5 7 10 20 25 30 40 50 55 58 60 70 80 100 200 300 400];%energies of proton projectile
  n = length(Ent);
  Zp = 1; %Z proton
  Za = 1; %Z hydrogen
  Io = 1; %13.6;%eV
  ao=5.3e-9; %cm

    function z = integrnd(w,Q) 
        F1=1./(1-exp(-2.*pi./sqrt(w-1)));
        F2=(Q+w/3)./Q;
        F3=(4*Q+(w-Q).^2).^3;
        FF1 = atan((2.*sqrt(w-1))./(2+Q-w));

    z = (F1.*F2.*exp((-2.*FF1)./sqrt(w-1)))./F3;

    end 
    for iii=1:n
    E=Ent(iii)
    vp=0.2*sqrt(E);
    Qmin = @(w) ((w)./(2.*vp)).^2;
    QSE = quad2d(@integrnd,1,2000,Qmin,2000,'MaxFunEvals',8000);

    saida=QSE %check the output

    cte =(2^9).*(Zp^2)./(vp.^2);%units of pi ao^2
    crosssection(iii) = cte.*QSE;
    end

    plot(Ent,crosssection,'*');
    end

Answer 1

[Mike Hosea]

I've only spent a little time looking at this. I get the warning but with passing the global error test. My impression is that the integral is just getting more difficult as E increases. There might be some precision difficulties with the way the integrand is coded. The expression for F2 seems problematic as Q-->0, which is happening as E is increased. More importantly, the expression for FF1 looks like it could swing wildly between positive and negative angles when Q is approximately equal to w-2. This would probably give QUAD2D fits. The latter problem can often be fixed by splitting the integral into two parts, each without discontinuity in the interior of the region of integration. However, in this case I wonder if that shoudn't be abs(2+Q-w) in the denominator. I don't know if that makes sense, but the integrator blasts right through then, and the resulting curve looks plausible.