From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: What am I doing wrong?
Date: Wed, 6 Apr 2011 19:36:05 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 26
Message-ID: <inifb5$49l$>
References: <> <inh86l$eof$>
Reply-To: <HIDDEN>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: 1302118565 4405 (6 Apr 2011 19:36:05 GMT)
NNTP-Posting-Date: Wed, 6 Apr 2011 19:36:05 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: comp.soft-sys.matlab:720464

"Roger Stafford" wrote in message <inh86l$eof$>...
> Paul Skoczylas <> wrote in message <>...
> > .........
> > 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)


 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