Advanced Mathematics and Mechanics Applications Using MATLAB phần 2 ppt

176: R2=[real(u2);imag(u2);z2];

177:

178: % Get curve properties from crvprp3d 179: [T,N,B,kappa]=crvprp3d(R1,R2);

180:

181: %==============================================

182:

183: function val=splined(xd,yd,x,if2) 184: %

185: % val=splined(xd,yd,x,if2)

186: % ~~~~~~~~~~~~~~~~~~~~~~~~ 187: %

188: % This function evaluates the first or second

189: % derivative of the piecewise cubic 190: % interpolation curve defined by the intrinsic

191: % function spline provided in MATLAB.If fewer 192: % than four data points are input, then simple

193: % polynomial interpolation is employed

194: % 195: % xd,yd - data vectors determining the spline

196: % curve produced by function spline

197: % x - vector of values where the first or 198: % the second derivative are desired

199: % if2 - a parameter which is input only if 200: % y’’(x) is required. Otherwise, y’(x)

201: % is returned.

202: % 203: % val - the first or second derivative values

204: % for the spline

205: % 206: % User m functions called: none

207:

208: n=length(xd); [b,c]=unmkpp(spline(xd,yd));

209: if n>3 % Use a cubic spline

210: if nargin==3, c=[3*c(:,1),2*c(:,2),c(:,3)]; 211: else, c=[6*c(:,1),2*c(:,2)]; end

212: val=ppval(mkpp(b,c),x);

213: else % Use a simple polynomial 214: c=polyder(polyfit(xd(:),yd(:),n-1));

215: if nargin==4, c=polyder(c); end 216: val=polyval(c,x);

217: end

218:

219: %=================================================

220:

© 2003 by CRC Press LLC

221: % function range=cubrange(xyz,ovrsiz)

222: % See Appendix B

2.8.2 Surface Properties

Surfaces are two-dimensional regions described parametrically as

R(u, v) = ˆıx(u, v) + ˆu(u, v) + kˆz(u, v)

where u and v are scalar parameters. This parametric form is helpful for generating a

grid of points on the surface as well as for computing surface tangents and the surface

normal. Holding v fixed while u varies generates a curve in the surface called a u

coordinate line. A tangent vector to the u-line is given by

gu = ∂R

∂u = ˆı

∂x

∂u + ˆ

∂y

∂u + kˆ ∂z

∂u.

Similarly, holding u fixed and varying v produces a v-line with tangent vector

gv = ∂R

∂v = ˆı

∂x

∂v + ˆ

∂y

∂v + kˆ ∂z

∂v .

Consider the following cross product.

gu × gv du dv = ˆn dS.

In this equation nˆ is the unit surface normal and dS is the area of a parallelogram

shaped surface element having sides defined by gu du and gv dv.

The intrinsic functions surf(X,Y,Z) and mesh(X,Y,Z) depict surfaces by showing

a grid network and related surface patches characterized when parameters u and v

are varied over constant limits. Thus, values

(uı, v) , 1 ≤ ı ≤ n , 1 ≤  ≤ m

lead to matrices

X = [x(uı, v)] , Y = [y(uı, v)] , Z = [z(uı, v)]

from which surface plots are obtained. Function surf colors the surface patches

whereas mesh colors the grid lines.

As a simple example, consider the ellipsoidal surface described parametrically as

x = a cos θ cos φ, y = b cos θ sin φ, z = c sin θ

where −π

2 ≤ θ ≤ π

2 , −π ≤ φ ≤ π. The surface equation evidently satisfies the

familiar equation

x

a

2

+

y

b

2

+

z

c

2

= 1

© 2003 by CRC Press LLC

for an ellipsoid. The function elipsoid(a,b,c) called with a = 2, b = 1.5, c = 1

produces the surface plot in Figure 2.18.

Many types of surfaces can be parameterized in a manner similar to the ellipsoid.

We will examine two more problems involving a torus and a conical frustum. Con￾sider a circle of radius b lying in the xz-plane with its center at [a,0,0]. Rotating

the circle about the z-axis produces a torus having the surface equation

x = [a + b cos θ] cos φ, y = [a + b cos θ] , sin φ, z = b sin φ

where −π ≤ θ ≤ π , −π ≤ φ ≤ π.

This type of equation is used below in an example involving several bodies. Let

us also produce a surface covering the ends and side of a conical frustum (a cone

with the top cut off). The frustum has base radius r b, top radius rt, and height h,

with the symmetry axis along the z-axis. The surface can be parameterized using an

azimuthal angle θ and an arc length parameter relating to the axial direction. The

lateral side length is

rs = h2 + (rb − rt)2 .

Let us take 0 ≤ s ≤ (rb + rs + rt) and describe the surface R(s, θ) by coordinate

functions

x = r(s) cos θ, y = r(s) sin θ, z = z(s)

where 0 ≤ θ ≤ 2π and

r(s) = s , 0 ≤ s ≤ rb

r(s) = rb + (rt − rb)(s − rb)

rs

, z = h(s − rb)

rs

, rb ≤ s ≤ (rb + rs)

r(s) = rb + rs + rt − r, z = h , (rb + rs) ≤ s ≤ (rb + rs + rt).

The function frus produces a grid of points on the surface in terms of r b, rt, h, the

number of increments on the base, the number of increments on the side, and the

number of increments on the top. Figure 2.16 shows the plot generated by frus.

An example called srfex employs the ideas just discussed and illustrates how

MATLAB represents several interesting surfaces. Points on the surface of an an￾nulus symmetric about the z-axis are created, and two more annuli are created by

interchanging axes. A pyramid with a square base is also created and the combina￾tion of four surfaces is plotted by finding a data range to include all points and then

plotting each surface in succession using the hold instruction (See Figure 2.16). Al￾though the rendering of surface intersections is not perfect, a useful description of a

fairly involved geometry results. Combined plotting of several intersecting surfaces

is implemented in a general purpose function surfmany. The default data case for

surfmany produces the six=legged geometry shown in Figure 2.17.

This section is concluded with a discussion of how a set of coordinate points can

be moved to a new position by translation and rotation of axes. Suppose a vector

r = ˆıx + ˆy + kˆz

© 2003 by CRC Press LLC

Spike and Intersecting Toruses

Figure 2.16: Spike and Intersecting Toruses

undergoes a coordinate change which moves the initial coordinate origin to (X o, Yo, Zo)

and moves the base vectors ˆı, ˆ, kˆ into eˆ1, eˆ2, eˆ3. Then the endpoint of r passes to

R = ˆıX + ˆY + kˆZ = Ro + ˆe1x + ˆe2y + ˆe3z

where

Ro = ˆıXo + ˆYo + kˆZo .

Let us specify the directions of the new base vectors by employing the columns of a

matrix V where we take

eˆ3 = V (:, 1)

norm[V (:, 1)] .

If V (:, 2) exists we take V (:, 1) × V (:, 2) and unitize this vector to produce eˆ2. The

triad is completed by taking eˆ1 = ˆe2 × eˆ3. In the event that V (:, 2) is not provided,

we use [1;0;0] and proceed as before. The functions rgdbodmo and rotatran

can be used to transform points in the manner described above.

© 2003 by CRC Press LLC

−5

0

5

−5

0

5

−8

−6

−4

−2

0

2

4

6

8

x axis

SEVERAL SURFACES COMBINED

y axis

z axis

Figure 2.17: Surface With Six Legs

© 2003 by CRC Press LLC

2

2.5

3

3.5

4

3

3.5

4

4.5

5

4

4.5

5

5.5

6

x axis

ROTATED AND TRANSLATED ELLIPSOID

y axis

z axis

Figure 2.18: Rotated and Translated Ellipsoid Surfaces

© 2003 by CRC Press LLC