Numerical integration of a 2d field with integration over parameters

17 ビュー (過去 30 日間)
dont panic
dont panic 2013 年 9 月 5 日
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 件のコメント
Roger Stafford
Roger Stafford 2013 年 9 月 5 日
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.)

サインインしてコメントする。

回答 (1 件)

Roger Stafford
Roger Stafford 2013 年 9 月 5 日
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 件のコメント
Roger Stafford
Roger Stafford 2013 年 9 月 5 日
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.

サインインしてコメントする。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by