Re: Surface integral
- To: mathgroup at smc.vnet.net
- Subject: [mg9854] Re: [mg9774] Surface integral
- From: Allan Hayes <hay at haystack.demon.co.uk>
- Date: Fri, 28 Nov 1997 05:36:10 -0500
- Sender: owner-wri-mathgroup at wolfram.com
Yves Verbandt wrote: > > Deat Mathematica experts, > > Can anyone help me with the calculation of a surface in 3D which is > given by a list of {x,y,z} values? Perhaps there is an adaptation of > ListIntegrate somewhere out there? > Yves: The area of a triangle with verticies a,b,c in any R^n is area[{a_,b_,c_}]:= Sqrt[Abs[Det[#.Transpose[#]]]]/2&[{b-a,c-a}] Take a matrix of points data=Table[{x,y,Sin[x-Sin[y]]}//N, {x,1,5},{y,1,6}]; Plot it ListPlot3D[Map[Last,data,{2}], AxesLabel-> {"x","y","z"}]; Cosider a point, p, in data (not along the bottom or the right of the matrix), with the neighboring points, e (east of p), s (south p} and se (south east of p). p--- e | | | | s--- se The sum of the areas of the triangles with vertex lists {p,s,e} and {se e,s} is given by the function area[{p_,s_,e_,se_}]:= area[{p,s,e}] + area[{se,s,e}] We shall replace each such p in the matrix with the quadruple {p,se,s,e} and then sum the values of area[{p,se,s,e}] p = Drop[#,-1]&/@Drop[data,-1]; s=Drop[#,-1]&/@Drop[RotateLeft[data],-1]; e=Drop[#,-1]&/@Drop[RotateLeft[data,{0,1}],-1]; se=Drop[#,-1]&/@Drop[RotateLeft[data,{1,1}],-1]; Make the matrix: data4 = MapThread[List,{p,s,e,se},2]; Make this into a list of quadruples {p,s,e,se} data4 = Flatten[data4,1] Now sum these areas Plus@@(area/@data4) 25.1413 Thes steps could be combined. We can compare this with the area of the related smooth surface: Clear[x,y] NIntegrate[ Evaluate[Sqrt[1+ D[Sin[x-Sin[y]],x]^2 + D[Sin[x-Sin[y]],y]^2]], {x,1,5},{y,1,6}] 25.6051 Another approximation is given by using the area of the parallelogram {p,s,e,s+e-p} instead of area[{p,s,e}] + area[{se,s,e}] This is the approximation underlying theintegral above. p = Drop[#,-1]&/@Drop[data,-1]; s=Drop[#,-1]&/@Drop[RotateLeft[data],-1]; e=Drop[#,-1]&/@Drop[RotateLeft[data,{0,1}],-1]; data4p= Flatten[MapThread[List,{p,s,e},2],1]; Plus@@(2 area/@data4p) 25.2283 -- Allan Hayes Mathematica Training and Consulting Leicester, UK hay at haystack.demon.co.uk http://www.haystack.demon.co.uk voice: +44 (0)116 271 4198 fax: +44 (0)116 271 4198