Integration using quadl with complex arguments

I am trying to integrate a function of the form:

quadl('s*real(R1)*besselj(0,s*rho),'0,Inf)

but I kept getting the following error:

??? Error using ==> inline.feval at 23 Not enough inputs to inline function.

Error in ==> quadl at 70 y = feval(f,x,varargin{:}); y = y(:).';

any clue on how I should sort this one out?

Your help would be greatly appreciated. I am putting my entire code below to put things in perspective:

%This program is an implementation of AKienle's 'Investigation of 2-layered %media with time-resolved reflectance Oct 1998' in the frequency domain.

%Elleesse

format long; ch='y'; global s;

% Enter inputs  
%  z = input('Enter z coordinate (mm) > ');
  rho=input('distance from the origin (mm) >');

% Specify constants
  b=rho;
  ua1=0.01;
  ua2=0.02;
  us1=1;
  us2=1;
  L0=10;
  f=110e6;
  c=0.299792458; %(mm/ps)
  v=c/1.4;
  om=2*pi*1e-12*f;
  D1=1./(3.*us1);
  D2=1./(3.*us2);
  dan12=1.4; %This is the refractive index mismatch between the scattering and non-scattering regimes

    if (dan12 > 1)
        A=504.332889-2641.00214*dan12+5923.699064*dan12.^2-7376.355814*dan12.^3+5507.53041*...
            dan12.^4-2463.357945*dan12.^5+610.956547*dan12.^6-64.8047*dan12.^7;
    elseif (dan12 <1)
        A=3.084635-6.531194*dan12+8.357854*dan12.^2-5.082751*dan12.^3+1.171382*dan12.^4;
    else
        A=1
    end

    % initialize value
    zb=2.*A.*D1;
    z0=1/us1;

      af1=sqrt(s^2+ua1./D1+1i.*om./(v*D1));
      af2=sqrt(s^2+ua2./D2+1i*om./(v*D2));

      % fi1 and R1 inherit their real and imaginary components from af1 and
      % af2. 
      fi1=(sinh(af1*(zb+z0))/(D1*af1)*(D1*af1*cosh(af1*L0)+D2*af2*sinh(af1*L0))/...
        (D1*af1*cosh(af1*(L0+zb))+D2*af2*sinh(af1*(L0+zb)))-sinh(af1*z0))/(D1*af1);

      %The derivative of fi1 is calculated and multiplied with D1. Fick's
      %Law was applied here (?)
      R1=-(sinh(af1*(zb+z0))*(D1*af1*sinh(af1*L0)+D2*af2*cosh(af1*L0))/...
         (D1*af1*cosh(af1*(L0+zb))+D2*af2*sinh(af1*(L0+zb)))) + cosh(af1*z0);

     rrlfi1 = real(fi1);
     imgnryfi1 = imag(fi1);
     rrlR1 = real(R1);
     imgnryR1 = imag(R1);

    w = quadl('s*rrlfi1*besselj(0,s*b)',0,10);
    x = quadl('s*imgnryfi1*besselj(0,s*b)',0,Inf);
    y = quadl('s*rrlR1*besselj(0,s*b)',0,Inf);
    a = quadl('s*imgnryR1*besselj(0,s*b)',0,Inf);

      % The integration is performed here using QUADL
        % sum_re=sum_re+besselj(0,s*rho)*re_R1(jj)*s;
        % sum_im=sum_im+besselj(0,s*rho)*im_R1(jj)*s;
        % sum_re1=sum_re1+besselj(0,s*rho)*re_fi1(jj)*s;
        % sum_im1=sum_im1+besselj(0,s*rho)*im_fi1(jj)*s;

 %---------------------------------------------------------------------
 % Revisions to the code

      %The formulae below gives the reflectance, R=R(rho)
      %sum_re=sum_re*si*0.306/(2*pi);
      %sum_im=sum_im*si*0.306/(2*pi);
      %sum_re1=sum_re1*0.118*si/(2*pi);
      %sum_im1=sum_im1*0.118*si/(2*pi);

    %  AC=sqrt((g+m)^2+(q+l)^2)
    %  phase=atan2(l+q,g+m)

    % You can't plot an AC vs s or phase vs s curve for that matter since
    % the program only gives one point by the time the AC and the phase are
    % calculated.

    AC=sqrt(((0.306*y)^2)+((0.306*a)^2))
    phase=-atan2((0.118*x)+(0.306*a),(0.306*y)+(0.118*w))
   ch = input('repeat the calculation? (y/n)', 's');

%----------------------------------------------------------------------

Answer 1

[Mike Hosea]

You cannot carry forward the dependence on 's' using simple variables. Let me restate Titus' suggestion in stronger words. Do not use a string input. Do not use an "inline" function. Instead, do use function handle syntax. I think this may look a little strange if you are not used to it, but in reality it is very easy. The beginning @(x) means that you are writing a function, and what follows is a going to be a function of x. Then just follow that with the expression you would evaluate on the command line in MATLAB. Any other variables that appear will have their current values substituted so that x will be the only unknown.

QUADL cannot integrate with infinite limits. You must use QUADGK.

So

w = quadl('s*rrlfi1*besselj(0,s*b)',0,10);

should be

w = quadgk(@(s)s.*rrlfi1(s).*besselj(0,s*b),0,10)

where previously we have defined

rrlfi1 = @(s)real(fil(s))

and before that defined

fil = @(s) (big long expression goes here, but it involves af1(s))

and before that we would have declared

af1 = @(s)sqrt(s.^2+ua1./D1+1i.*om./(v*D1));

Now, I'm assuming all thos other variables were scalar constants known at the time af1 is defined. Make sure where s is involved, you use .*, ./, or .^ instead of *, /, or ^ because your function has to work with vector input for s for QUADGK (or QUADL, for that matter) to work on it.

[Elleesse]

Thank you so much for responding, Mike. I ran the code again using quadgk and it worked handsomely for a considerable range of rho.

Answer 2

[Titus Edelhofer]

Hi,

you would need to pass the parameters to the inline functions. I'd suggest to use anonymous functions instead, i.e., replacing

w = quadl('s*rrlfi1*besselj(0,s*b)',0,10);

by

w = quadl(@(s) s*rrlfi1*besselj(0,s*b), 0, 10);

BTW: When I tried the code, rrlfi1 was empty, so the quadl would fail here anyway. But I didn't look why rrlfi1 was empty ...

Titus

[Elleesse]

Thank you so much for responding. I really appreciate it.

You're right: rrlfi1 was empty and I couldn't figure it out either. Let me describe what I am trying to do:

* These are the functions I want to integrate with respect to the variable 's' using quadl (w,x,y and a above) from 0 to infinity:

s*real(R1)*besselj(0,s*rho)

s*real(fi1)*besselj(0,s*rho)

s*imag(R1)*besselj(0,s*rho)

s*imag(fi1)*besselj(0,s*rho)

*R1 and fi1 have an implicit dependence on the variable 's' by their dependence on af1 and af2.

*Being new to MATLAB, I thought that by virtue of that dependence, and the fact that I had declared 's' to be a global variable, I can just pass R1 and fi1 to the functions I want to integrate and quadl would take care of the integration.

*But, as you yourself found out, it just didn't work out. Alternatively, if I do assign 's' at the beginning of the code with a numerical value, then the integration wouldn't integrate the function with respect to 's' or with respect to any variable for that matter.

For rho = 8.00063, AC and phase should equal to:

AC = 3.41E-04

phase = 0.15

or somewhat close to that value.

Your insights would be greatly valued.