I am getting the following error message and would like to know how I can remedy the integral error. Thanks.

Error using integralParseArgs (line 11)
Expected a string scalar or character vector for the parameter name.

Error in integral (line 87)
opstruct = integralParseArgs(varargin{:});

Error in Ass_2_AcousticOceans (line 120)
                EJP(iy,ix,im)=integral(@(n2)f,n2,-inf,inf);

clear all
close all
clc

c1=1500;
c2=1800;
row1=1000;
row2=1200;
f=150;
w=2*pi*f;
h=50;
r=500;
k1=w/c1;
k2=w/c2;
b12=row1/row2;
Q=1;
m=5;
z=20;
zprime=40;
alphac=acos(c1/c2);

%Smz1=sin(n1*z);                   %mode shape at fluid channel
%Smz2=sin(n1*h)exp(-i*n2m(z-h));    %mode shape at fluid basement

if z<h                              %set z greater and lesser
    zl=z;
    zg=zprime;
else
    zl=zprime;
    zg=z;
end

% U=sqrt(1-(m*pi)^2/(k1*he*sin(alphac))^2);
% S1=sin(m*pi*z/he);                                  %propagating mode shape  0<z<h
% S2=sin(m*pi*z/he).*exp(-k1*(z-h)*sin(alphac)*U);    %evanescent mode shape   z>h

phi=zeros(h,r,m);
EJP=zeros(h,r,m);
Res=zeros(h,r,m);

for im=1:1:m
    for ix=1:1:r
        for iy=1:1:h
            he=h*(1+1/(b12*k1*h*sin(alphac)));
            n1=(im*pi)/he;
            n2=i*k1*sin(alphac)*(1-(im^2*pi^2)/(h^2*he^2*sin(alphac)))^0.5;
            p=-k1*(1-(im^2*pi^2)/(k1^2*he^2))^0.5;
            D=n1*cos(n1*h)+i*b12*n2*sin(n1*h);
            F1=(sin(n1*zl)/n1)*(n1*cos(n1*(h-zg))+i*b12*n2*sin(n1*(h-zg)))/D;
            F2=(sin(n1*zl)*exp(-i*n2*(zg-h)))/D;
            
            if h<=h
                Res(iy,ix,im)=(n1*sin(n1*z)*sin(n1*zprime))/...
                    (p*(n1*h-sin(n1*h)*cos(n1*h)-b12^2*sin(n1*h)*tan(n1*h)))*exp(i*p*ix);
                f=n2/(k2^2-n2^2)^0.5*F1*exp(-i*(k2^2-n2^2)^0.5*ix);
                EJP(iy,ix,im)=integral(@(n2)f,n2,-inf,inf);
                else 
                Res(iy,ix,im)=b12*(n1*sin(n1*h)*sin(n1*zprime)*exp(-i*n2(z-h)))/...
                    (p(n1*h-sin(n1*h)*cos(n1*h)-b12^2*sin(n1*h)*tan(n1*h)))*exp(i*p*ix);
                f=n2/(k2^2-n2^2)^0.5*F2*exp(-i*(k2^2-n2^2)^0.5*ix);
                EJP(iy,ix,im)=b12*integral(@(n2)f,n2,-inf,inf);
            end
            phi(iy,ix,im)=(Q/(2*pi*ix^0.5))*(2*pi*i*Res(iy,ix,im)+EJP(iy,ix,im));
        end
    end
end

phi=sum(phi,3);
pcolor(h,x,abs(psi));

[KSSV]

YOu need to define a anonymous function using @..there is no such function???? Read about integral function carefully.

[Brandon Johnson]

I define it within the integral function. That should work as well. If i define it outside the integral then it gives the same error message.

Answer 1

[Walter Roberson]

EJP(iy,ix,im)=integral(@(n2)f,n2,-inf,inf);

asks to run integral with anonymous function @(n2)f . The anonymous function takes a single parameter n2 and ignores it and returns whatever f is, which is an error because you need to return something the same size as the input.

The second parameter is the lower bound, and is given as n2 . The n2 that would be referred to would be the one constructed at

n2=i*k1*sin(alphac)*(1-(im^2*pi^2)/(h^2*he^2*sin(alphac)))^0.5;

This would not be the same n2 as the n2 referred to inside @(n2)f where it is used as a parameter name shadowing the assigned-to n2.

The third parameter is the upper bound, and is given as -inf . An upper bound of -inf is unusual, but not prohibited. As -inf is going to be less than whatever n2 was calculated at, you would get the negative from what you would have calculated if you had used -inf as the lower bound and n2 as the upper bound.

The fourth parameter, when present, must be a keyword for a keyword parameter pair. Instead you give +inf.

[Brandon Johnson]

I've redefined the integral like this and that worked better as I want to valuate the integral from -inf to inf, with relation to the variable n2. Now I'm getting a whole new set of errors and am not to sure what to make of them,

EJP(iy,ix,im)=integral(@(n2)f,-inf,inf);

EJP(iy,ix,im)=b12*integral(@(n2)f,-inf,inf);

Error using integralCalc/finalInputChecks (line 526)
Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the 'ArrayValued'
option to true.

Error in integralCalc/iterateScalarValued (line 315)
                finalInputChecks(x,fx);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 103)
        [q,errbnd] = vadapt(@minusInfToInfInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

Error in Ass_2_AcousticOceans (line 56)
                EJP(iy,ix,im)=integral(@(n2)f,-inf,inf);
 

[Walter Roberson]

EJP(iy,ix,im)=integral(@(n2)f,-inf,inf);

asks to run integral with anonymous function @(n2)f . The anonymous function takes a single parameter n2 and ignores it and returns whatever f is, which is an error because you need to return something the same size as the input.

Your f is defined as an expression

f=n2/(k2^2-n2^2)^0.5*F1*exp(-i*(k2^2-n2^2)^0.5*ix);

Note that the n2 there will have nothing to do with the n2 of @(n2)f . You had calculated

n2=i*k1*sin(alphac)*(1-(im^2*pi^2)/(h^2*he^2*sin(alphac)))^0.5;

before that point, giving a numeric value, and that numeric value would have been used to calculate a numeric f, and the value n2 would have lost its identity inside f. You are not building a formula named f in terms of n2, or a function named f in terms of n2: you calculated a numeric value, just like

a = 1;
b = 2;
c = a + b;
a = 5;

then c does not suddenly become 5 + 2 = 7: it stays as 3, since the a+b is not creating a formula or a function.

Now, perhaps you want

f = @(n2) n2 ./ (k2.^2-n2.^2).^0.5 .* F1 .* exp(-i.*(k2.^2-n2.^2).^0.5 .* ix);

but if you think you want that, remember the numeric value F1 was build from the numeric n2 that you assigned earlier.

Is it plausible that you want to go back all those levels and substitute in the n2 from the integral(), with n2 ranging from -inf to +inf ? No, it is not plausible. You define n2 in terms of i times something, and since you have not assigned to i, i stands for the imaginary unit, sqrt(-1), so unless the (1-(im^2*pi^2)/(h^2*he^2*sin(alphac))) is certain to non-positive so that the sqrt() of it is going to generate a sqrt(-1) that you are canceling out with the i* then we must assume that you want n2 to be complex valued, which would disagree with sweeping n2 over the real range -inf to +inf.

.... Basically it looks to us observers as if your formulas are messed up for what you are asking to integrate.

[Brandon Johnson]

I've fixed it.. The formulas are correct I just needed to set the n2 for the Res function and it being intergates from -inf to inf was correct for that function. I also needed to extend the Relerror for the integration.

clear all
close all
clc

c1=1500;
c2=1800;
row1=1000;
row2=1800;
f=150;
w=2*pi*f;
h=50;
r=50;
k1=w/c1;
k2=w/c2;
b12=row1/row2;
Q=1;
m=5;
z=20;
zprime=40;
alphac=acos(c1/c2);

%Smz1=sin(n1*z);                   %mode shape at fluid channel
%Smz2=sin(n1*h)exp(-i*n2m(z-h));    %mode shape at fluid basement

if z<h                              %set z greater and lesser
    zl=z;
    zg=zprime;
else
    zl=zprime;
    zg=z;
end

% U=sqrt(1-(m*pi)^2/(k1*he*sin(alphac))^2);
% S1=sin(m*pi*z/he);                                  %propagating mode shape  0<z<h
% S2=sin(m*pi*z/he).*exp(-k1*(z-h)*sin(alphac)*U);    %evanescent mode shape   z>h

phi=zeros(h,r,m);
EJP=zeros(h,r,m);
Res=zeros(h,r,m);

for im=1:1:m
    for ix=1:1:r
        for iy=1:1:h
            he=h*(1+1/(b12*k1*h*sin(alphac)));
            n1=(im*pi)/he;
            n2=i*k1*sin(alphac)*(1-(im^2*pi^2)/(h^2*he^2*sin(alphac)))^0.5;
            p=-k1*(1-(im^2*pi^2)/(k1^2*he^2))^0.5;
            D=n1*cos(n1*h)+i*b12*n2*sin(n1*h);
            F1=(sin(n1*zl)/n1)*(n1*cos(n1*(h-zg))+i*b12*n2*sin(n1*(h-zg)))/D;
            F2=(sin(n1*zl)*exp(-i*n2*(zg-h)))/D;
            
            if h<=h
                Res(iy,ix,im)=(n1*sin(n1*z)*sin(n1*zprime))/...
                    (p*(n1*h-sin(n1*h)*cos(n1*h)-b12^2*sin(n1*h)*tan(n1*h)))*exp(i*p*ix);
                f=n2/(k2^2-n2^2)^0.5*F1*exp(-i*(k2^2-n2^2)^0.5*ix);
                EJP(iy,ix,im)=integral(@(n2)f,-inf,inf,'ArrayValued',true,'AbsTol',1e-3,'RelTol',1);
                else 
                Res(iy,ix,im)=b12*(n1*sin(n1*h)*sin(n1*zprime)*exp(-i*n2(z-h)))/...
                    (p(n1*h-sin(n1*h)*cos(n1*h)-b12^2*sin(n1*h)*tan(n1*h)))*exp(i*p*ix);
                f=n2/(k2^2-n2^2)^0.5*F2*exp(-i*(k2^2-n2^2)^0.5*ix);
                EJP(iy,ix,im)=b12*integral(@(n2)f,-inf,inf,'ArrayValued',true,'AbsTol',1e-3,'RelTol',1);
            end
            phi(iy,ix,im)=(Q/(2*pi*(ix^0.5+abs(iy-z)^2)^0.5))*(2*pi*i*Res(iy,ix,im)+EJP(iy,ix,im));
        end
    end
end

phi=sum(phi,3);
phi =abs(reshape(phi,[],size(phi,2),1));
PP = 20*log10(abs(phi));

figure (1)
imagesc(1:1:r,1:1:h,PP);hold on
title('Pressure field (dB) ','fontweight','b','FontSize',16);
xlabel('Range (meter)','fontweight','b','FontSize',16);
ylabel('Depth (meter)','fontweight','b','FontSize',16);
set(gca,'Ydir','Reverse','fontweight','b','FontSize',16);%set(gca,'Ydir','reverse')

hold off
colormap jet
colorbar
caxis([30 100])

figure (2)
pcolor(1:1:r,1:1:h,abs(phi));
set(gca,'Ydir','reverse','fontweight','b','FontSize',10);
shading interp; 
colormap(jet); 
colorbar;
title('Particle Velocity Potential Field (1/m^3)');