Numerical integration of a 2d field with integration over parameters

Hello,

i have tried to numerically integrate a funktion f(r,phi,z,R,Z) of which R, and Z are the field variables of my 2d field and r,phi and z are parameters which i want to integrate numerically.

My code looks like this :

%Function

a=0.0025; L=0.05; phis=0;

f=@(r,phi,z,R,Z) (sin(phis-phi).*r)/(R.^2-2.*r.*R.*cos(phis-phi)+r.^2+(Z-z).^2).^(3/2);

%integration boundaries

r1=0; r2=a; phi1=0; phi2=2*pi; z1=-L/2; z2=L/2;

[R,Z] = meshgrid(-0.05:s:0.05); %field variables mp=zeros(size( R ));

for k = 1:numel( R )

    g = @(r,phi,z) f(r,phi,z,R(k),Z(k));
    mp(k) = triplequad(g,r1,r2,phi1,phi2,z1,z2);

end

surf(R,Z,mp)

But I always get a error message.

Which looks like this:

??? Error using ==> quad at 79 The integrand function must return an output vector of the same length as the input vector.

Error in ==> triplequad>innerintegral at 63 Q(i) = quadf(intfcn, xmin, xmax, tol, trace, y(i), z, varargin{:});

Error in ==> dblquad>innerintegral at 73 fcl = intfcn(xmin, y(1), varargin{:}); %evaluate only to get the class below

Error in ==> quad at 76 y = f(x, varargin{:});

Error in ==> dblquad at 53 Q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ...

Error in ==> triplequad at 47 Q = dblquad(@innerintegral, ymin, ymax, zmin, zmax, tol, quadf, intfcn, ...

Error in ==> pot_keep_parameter at 41 mp(k) = triplequad(g,r1,r2,phi1,phi2,z1,z2);

Thank you for your help

[Roger Stafford]

The 'triplequad' function may be complaining about the absence of a dot in the division inside f. It should presumably be:

 f=@(r,phi,z,R,Z) (sin(phis-phi).*r)./(R.^2-2.*r.*R.*cos(phis-phi)+r.^2+(Z-z).^2).^(3/2);

Note the added 'dot' in front of the '/' symbol. This will assure that the vector output will be the same length as the vector input, since the first variable, r, may be given by 'triplequad' in the form of a vector, not a scalar. (We wouldn't want matrix division to occur at this point in any case.)

Answer 1

[Roger Stafford]

If your triple integral is intended to be the volume integral over the defined limits of cylindrical coordinates, r, phi, and z, you will need to include the Jacobian factor, r, in your integrand to make it valid. I don't see any evidence of it there. The 'r' that is there is apparently needed to obtain the force component orthogonal to the R,Z plane and that is in addition to the Jacobian factor.

Of course this doesn't account for your error message. In the text of your for-loop I get a peculiar symbol after the 'numel' instead of numel(mp). It looks like a registered trademark symbol. What is that? Can you please show us the entire error message?

[dont panic]

Thank you for your answer but i have included the r i : (sin(phis-phi).*_r_)

you are right is it supposed to be a volume integral in cylindrical coordinates and I#m integrating a potential of a pointsource over the volume of a body and want to calculate the potential field. The parameters for the source are:r,phi and z and that are the ones i integrate and the coordinates for the field are R and Z while phis is set to be 0 due to axial symetry.

[Roger Stafford]

If I am not mistaken, that r you have in the formula for f is needed for computing the cosine of the angle between the line connecting the (R,Z) point with the (r,phi,z) point and the normal to the plane of (R,Z) which would be:

 cosine of angle =
 (sin(phis-phi).*r)/(R.^2-2.*r.*R.*cos(phis-phi)+r.^2+(Z-z).^2).^(1/2)

and the coulomb law inverse distance squared factor would be the rest:

 1./(R.^2-2.*r.*R.*cos(phis-phi)+r.^2+(Z-z).^2).^1

The Jacobian factor is needed in addition to this to make the differential correction:

 dV = dx*dy*dz = r*dphi*dr*dz

where dV refers to volume differential. Otherwise you would receive the wrong answer for the volume integral.

[Roger Stafford]

I would like to add that that does not look like an electrical potential you are integrating. With a potential I would expect to see

 1/(R.^2-2.*r.*R.*cos(phis-phi)+r.^2+(Z-z).^2).^(1/2)

and not the 3/2 power, and also not accompanied by the factor

 sin(phis-phi).

What you have given in f is what I would expect to see with an integration of the component of electrical force normal to the R,Z mesh plane of a unit density charge distribution throughout the volume as exerted at a point in the mesh.