Numerical integration of a 2d field with integration over parameters
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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
1 Kommentar
Roger Stafford
am 5 Sep. 2013
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.)
Antworten (1)
Roger Stafford
am 5 Sep. 2013
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?
3 Kommentare
Roger Stafford
am 5 Sep. 2013
Bearbeitet: Roger Stafford
am 5 Sep. 2013
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
am 5 Sep. 2013
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.
Siehe auch
Kategorien
Mehr zu Surface and Mesh Plots finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!