Thread Subject:

What am I doing wrong?

Please forgive the horribly meaningless title... I couldn't think of a more descriptive title that would be meaningful.

I am trying replicate some work in an old paper. The paper is by Jaeger and Chamalaun and is online at http://adsabs.harvard.edu/full/1966AuJPh..19..475J I am trying to replicate the integral in Equation 10, which refers to Eq 11. (Click on the link for p 477.)

For the integrand, I have defined the functions listed below. I believe there's nothing wrong with my integration routines, so for now, I'm just asking if I've misinterpreted the paper for the description of the integrand (or some similar mistake). The odd thing is that it seems to work pretty well for some values of tau and B, but is wildly out to lunch for others.

Note that I have full confidence in the values published in the paper. I ran a numerical inverse Laplace transform on them and every single one matched exactly what was in the paper when rounded to the same number of decimals. But I'd like to be able to do a numerical integration at least as well as they could do in 1966!

Here are the functions--the first one is the one I'm calling from my integration routine. It's set up for u to be a vector, but tau and B are scalars. It calls the others.

function I=calcphiI1(u,tau,B);
   I=exp(-tau*u.^2).*(tau./(u.^2+B^2)+1./(u.^2+B^2).^2).*u.*(pi/2+arg((u.*Hankel1(1,u)+B*Hankel1(0,u))/B));

function H=Hankel1(a,x)
H=besselj(a,x)+i*bessely(a,x);

function a=arg(z);
x=real(z);
y=imag(z);

if x>0 | y~=0;
    a=2*atan(y./(sqrt(x.^2+y.^2)+x));
elseif x<0 & y==0;
    a=pi;
elseif x==0 &y==0
    a=NaN;
else;
    a=NaN;
end;


%%%

Thanks for any help anyone can provide!

-Paul

On Apr 6, 5:22 am, Paul Skoczylas <pa...@cfertech.com> wrote:
> Please forgive the horribly meaningless title...  I couldn't think of a more descriptive title that would be meaningful.
>
> I am trying replicate some work in an old paper.  
...
> Note that I have full confidence in the values published in the paper.  I ran a numerical inverse Laplace transform on them and every single one matched exactly what was in the paper when rounded to the same number of decimals.  But I'd like to be able to do a numerical integration at least as well as they could do in 1966!

Forget about that! The guys who worked with these kinds of
things in 1966 were those who first of all had the talents
needed to work with engineering, the motivation to work *hard*
to get that kind of job, and also had the tutors to match.

The fact that you ask this question in a matlab newsgroup,
and not in a numerical maths newsgrouop, means that you don't
have what it takes to come anywhere near succeeding with such
kinds of things.

Rune

On Apr 6, 4:41 pm, Rune Allnor <all...@tele.ntnu.no> wrote:
> On Apr 6, 5:22 am, Paul Skoczylas <pa...@cfertech.com> wrote:
>
> > Please forgive the horribly meaningless title...  I couldn't think of a more descriptive title that would be meaningful.
>
> > I am trying replicate some work in an old paper.  
> ...
> > Note that I have full confidence in the values published in the paper.  I ran a numerical inverse Laplace transform on them and every single one matched exactly what was in the paper when rounded to the same number of decimals.  But I'd like to be able to do a numerical integration at least as well as they could do in 1966!
>
> Forget about that! The guys who worked with these kinds of
> things in 1966 were those who first of all had the talents
> needed to work with engineering, the motivation to work *hard*
> to get that kind of job, and also had the tutors to match.
>
> The fact that you ask this question in a matlab newsgroup,
> and not in a numerical maths newsgrouop, means that you don't
> have what it takes to come anywhere near succeeding with such
> kinds of things.
>
> Rune

Rune

Back in those days, they probably also had Fortran, so they had a
language that was designed for such problems. ;-)

Paul Skoczylas <pauls@cfertech.com> wrote in message <ee902dac-f320-44bf-ae32-e6d5cabf7713@glegroupsg2000goo.googlegroups.com>...
> .........
> For the integrand, I have defined the functions listed below. I believe there's nothing wrong with my integration routines, so for now, I'm just asking if I've misinterpreted the paper for the description of the integrand (or some similar mistake). The odd thing is that it seems to work pretty well for some values of tau and B, but is wildly out to lunch for others.
> ..........
> function H=Hankel1(a,x)
> H=besselj(a,x)+i*bessely(a,x);
>
> function a=arg(z);
> x=real(z);
> y=imag(z);
>
> if x>0 | y~=0;
> a=2*atan(y./(sqrt(x.^2+y.^2)+x));
> elseif x<0 & y==0;
> a=pi;
> elseif x==0 &y==0
> a=NaN;
> else;
> a=NaN;
> end;
> .........
- - - - - - - - - -
  I will hazard the guess that your trouble lies in the meaning of the 'arg' function. As u advances from 0, the phase angle of J0(u)+i*Y0(u) begins at -pi/2 and moves through the fourth quadrant into the first, then the second, the third, back to the fourth, on to the first, and so forth, always winding around in a clockwise direction of increasing phase. Something similar occurs with u/B*(J1(u)+i*Y1(u)). If continuity is to be preserved, the phase of the sum of these must be a quantity that does not take the backward jumps that would keep it within the bounds -pi to pi, and therefore you cannot calculate it as you have done in your 'arg'. The statement in the paper that the psi(u) ~ u as u approaches infinity is a clue that this is actually the case.

  You would need to define it in a very different manner perhaps using something like the 'unwrap' function to assist in this, though I am not quite sure offhand how you would do it in conjunction with your quadrature procedure.

  (By the way, Matlab's 'atan2' function does just what your currently defined 'arg' does. However, as I say, the 'arg' you actually need would presumably have to be different from this.)

Roger Stafford

"Roger Stafford" wrote in message <inh86l$eof$1@fred.mathworks.com>...
> Paul Skoczylas <pauls@cfertech.com> wrote in message <ee902dac-f320-44bf-ae32-e6d5cabf7713@glegroupsg2000goo.googlegroups.com>...
> > .........
> > For the integrand, I have defined the functions listed below. I believe .....
> I will hazard the guess that your trouble lies in the meaning of the 'arg' .....
- - - - - - - - - -
  Thinking about your problem again, I can see no other practical way of computing the psi(u) function than selecting a quadrature method that uses fixed, closely-spaced points of u, rather than attempting to define it as a general function as is called for with 'quad' and other routines like it. That way you can use the Matlab 'unwrap' function to avoid unwanted jumps of -2*pi. Letting u be a vector of u-values starting at (or awfully close to) zero and going to large enough values to approximate infinity for integration purposes, you can compute psi this way:

 psi = unwrap(atan2(u/B.*J1(u)+J0(u),-u/B.*Y1(u)-Y0(u)));

These would all satisfy the equations

 M(u)*cos(psi(u)) = -u/B.*Y1(u)-Y0(u)
 M(u)*sin(psi(u)) = u/B.*J1(u)+J0(u)

where

 M(u) = sqrt((u/B.*Y1(u)+Y0(u))^2+(u/B.*J1(u)+J0(u))^2).

Psi(u) should start at zero and continue increasing indefinitely without any discontinuous jumps. Of course you must make sure that the spacing in u is always close enough that 'unwrap' would not choose the wrong multiple of 2*pi, (but in fact they should be much closer than that to ensure integration accuracy.)

  For quadrature you can either use very closely spaced points with Matlab's 'trapz' or else use one of the higher order FEX routines with fewer points. Since the Bessel functions are smooth, one of the latter would probably work very well with, say, a third or fourth order routine.

  The reason for avoiding discontinuities in psi is of course that formula (10) was obtained using integration by parts and the validity of integration by parts depends very much on the continuity of psi and its derivative.

Roger Stafford