Numerical Optimization: Finding Local Extreme Values using Hunt Algorithm, Study notes of Mathematics

Download Numerical Optimization: Finding Local Extreme Values using Hunt Algorithm and more Study notes Mathematics in PDF only on Docsity! Some Calculus III Notes Carl Eberhart Intersection of 3 cylinders –1–0.500.51 x –1 –0.5 0 0.5 1 y –1 –0.5 0 0.5 1 z January 5, 2009 Extreme values > with(plots): #execute this line first 1 least squares fit If you have 4 points in space ( xi, yi, zi), i = 1..4, there probably isn’t a plane which contains all of them. However, you can find the plane z = ax + by + c which best fits the point in the sense that g(a, b, c) = ∑4 i=1 (axi + b yi + c − zi)2 is as small as possible. There is only one critical point for this function It is the solution of the system 3 linear equations in a,b and c ga = 0, gb = 0, gc = 0 > data:=([[0,0,2],[2,0,4],[2,3,8],[0,3,5]]): > matrix(data);      0 0 2 2 0 4 2 3 8 0 3 5      make up the function g > g := (c - 2)^2 + (2*a+c -4)^2 + (2*a+3*b+c-8)^2 + (3*b+c-5)^2; g := (c − 2)2 + (2 a + c − 4)2 + (2 a + 3 b + c − 8)2 + (3 b + c − 5)2 find its critical point. Since g is a sum of squares, it will have a minimum at the critical point. > sol:=solve({diff(g,a),diff(g,b),diff(g,c)},{a,b,c}); sol := {a = 5 4 , b = 7 6 , c = 7 4 } If we put these values for a, b, and c into g and then compute sqrt(g(a,b,c)/4) we get a measure of the goodness of fit. > goodnessoffit:=sqrt(subs(sol,g)/4); goodnessoffit := 1 4 this says on average the measured value is within 1/4 of the ’correct’ value. we can draw the points and the plane to see how well they fit. > display (plot3d(subs(sol,a*x+b*y+c),x=-1..3,y=-1..4,color=grey,style=wireframe ),pointplot3d({op(data)},symbol=cross,color=red),axes=boxed); 1 –0.5 0 0.5 1 1.5 x –0.5 0 0.5 1 1.5 y 6 7 8 9 10 Play around with hunt and showhunt. > plot3d(4*x*y+1-x^4-y^4,x=-2..2,y=-2..2,view=[-2..2,-2..2,-1..6]); > showhunt(4*x*y+1-x^4/2-y^4,.1,.1,.2,10,.03,4); –2 0 2 4 x –3 –2 –1 0 1 2 3 4 y –300 –200 –100 0 3 finding critical points algebraically and applying the 2nd deriva- tive test. 4 Here is a function to study. > f:=x^4+y^4-4*x*y + 1; f := x4 + y4 − 4x y + 1 Calculating the first partial derivatives of f. > fx := diff(f,x); fy:=diff(f,y); fx := 4x3 − 4 y fy := 4 y3 − 4x Using solve to find the critical points. We see that there are only three. The other two are complex pairs > evalf(solve({fx,fy},{x,y})); {x = 0., y = 0.}, {x = 1., y = 1.}, {x = −1., y = −1.}, {x = 1. I, y = −1. I}, {x = 0.7071067812 + 0.7071067812 I, y = −0.7071067812 + 0.7071067812 I} You can also use fsolve to find the critical points. > s1:=fsolve({fx,fy},{x,y},{x=-2..2,y=-2..2}); s1 := {x = 0., y = 0.} > s2:=fsolve({fx,fy},{x,y},{x=-2..2,y=-2..2},avoid={s1}); s2 := {x = −1.000000000, y = −1.000000000} > s3:=fsolve({fx,fy},{x,y},{x=-2..2,y=-2..2},avoid={s1,s2}); s3 := {x = 1.000000000, y = 1.000000000} Use the second derivative test on the 3 critical points > fxx:=diff(fx,x); fyy:=diff(fy,y); fxy:=diff(fx,y); fxx := 12x2 fyy := 12 y2 fxy := −4 > discrcim:=fxx*fyy-fxy^2; discrcim := 144x2 y2 − 16 We see that at (0,0) the discriminant is <0 so there is no local extreme. At (1,1) and (-1,-1) however, the discriminant is >0 and fxx >0 so the function has local maxima at (1,1) and (-1,-1). 4 maximizing the volume of a box with a half top, and a quarter front A box is to be made from 10 square feet of metal so that its top is only half there and its front is 3/4 missing. Find the maximum volume of the box. 5 Maximimize V = xyz subject to 10 = 3/2 xy + 2yz + 7/4 xz first solve the constraint for one of the variables and reduce V to a function of 2 variables. > constraint:=10 = 3/2* x*y + 2*y*z + 7/4 *x*z; constraint := 10 = 3 2 x y + 2 y z + 7 4 x z > sol := solve(constraint,z); sol := −2 (−20 + 3x y) 8 y + 7x > V := x*y*sol; V := −2x y (−20 + 3x y) 8 y + 7x Find critical points > Vx := diff(V,x); Vy:=diff(V,y); Vx := −2 y (−20 + 3x y) 8 y + 7x − 6x y 2 8 y + 7x + 14x y (−20 + 3x y) (8 y + 7x)2 Vy := −2x (−20 + 3x y) 8 y + 7x − 6x 2 y 8 y + 7x + 16x y (−20 + 3x y) (8 y + 7x)2 These are ugly. Lets simplify > Vx:=simplify(Vx); Vx := −2 y 2 (−160 + 48x y + 21x2) (8 y + 7x)2 > Vy:=simplify(Vy); Vy := −4x 2 (−70 + 12 y2 + 21x y) (8 y + 7x)2 By inspection, we see that Vx and Vy will = 0 when x = 0 = y. However these values are not in the domain of V So we next look for the zero’s of the other two factors of the top of Vx and Vy > eqns:={-160+48*x*y+21*x^2=0,-70+12*y^2+21*x*y=0}; eqns := {−160 + 48x y + 21x2 = 0, −70 + 12 y2 + 21x y = 0} We can draw these two curves in the plane using implicitplot. > implicitplot(eqns,x=0..10,y=0..10); 6 there isn’t a quadratic example of this. But if we take example and add xˆ3 to it, this turns (0,0) into a ’hanging valley’ and does not alter the value of D. > display( pointplot3d([0,0,0],symbol=circle,color=blue), spacecurve([t,t,t^3],t=-1..1, color=red), > plot3d(x^2-2*x*y+y^2+x^3,x=-1..1,y=-1..1),axes=boxed,orientation =[73,57]); –1–0.500.51 x –1 –0.5 0 0.5 1 y 0 1 2 3 4 5 6 finding extreme values of a continuous function on a closed and bounded set in the plane Here is an ’algorithm’ for finding the extreme values of a continuous function f on a closed and bounded domain D. 1. Identify the boundary of D and the interior of D. The interior of D is the set of points (a,b) in D for which there is a positive ε so that if (x,y) is within ε of (a,b), then (x,y) lies in D. The boundary of D is the points of D that are not in the interior of D. Those points have the property that there are points arbirarily close to them which are not in D. We assume that f is differentiable on the interior of D. 2. Find the critical points of f in the interior of D, and tabulate the function there. 3. Find the extreme values of f on the boundary of D. The boundary of D is usually a curve of some sort. So one way to optimize f on the boundary would be to parameterize the boundary with a function r(t) and then optimize the composition f r. Note 1: r might be piecewise defined. If the boundary is complicated, you might have to break it into pieces and parameterize the pieces. Note 2: Sometimes it is useful to draw a few gradient vectors for f to see how f is increasing. This might help narrow down the search for maximum and minimum values. So if all the gradient vectors are pointing into the interior away from boundary at a certain place, probably there is not a maximum value of f on the boundary there, because the gradient vectors point in the direction of (greatest) increasing f values. 9 Examples: 1 Optimize (ie find all extreme values of) f(x,y) = x y + x2 on D = the triangular region with vertices (0,2), (0,-2), and (1,-2) The boundary of D consists of the three segments connecting the vertices and the interior is all points inside the triangle. One checks that fx = 0 = fy only at (0,0) which is not in the interior of D. Thus the maximum value of f on D must lie on the boundary of D. We can find the extreme values on each edge separately: on the edge (0,2) to (0,-2) f(x,y) = 0 on the edge (0,2) to (1,-2), y = 2 − 4x so we want to optimize f(x, y) = f(x, 2 − 4x) = x (2 − 4x) + x2 = 2x − for 0 ≤ x ≤ 1 the derivative 2 − 6x = 0 when x = 1/3, so the extreme values on this interval must be among the values f(0,2) = 0, f(1/3,2/3) = 1/3 and f(1,-2) = -1. on the edge (0,-2) to (1,-2), f(x,y) = f(x,-2) = −2x + x2 for x between 0 and 1. the derivative -2 + 2x = 0 when x = 1, so the extreme values on this interval are f(0,-2) = 0 and f(1,-2) = -1 Combining these efforts, we see that f(1/3,2/3) = 1/3 is the maximum value and f(1,-2) = -1 is the minimum value of f(x,y) on the given triangular region. We can check our work visually in Maple > with(plots): This is a domain restriction function. It returns (x,y) if (x,y) is in the domain of f, otherwise it returns (0,0) (so that f(r(x,y)) = f(x,y) if (x,y) is in the domain of f and f(r(x,y)) = 0 otherwise. > r:= proc(x,y) if x>= 0 and y >=-2 and y <= 2 - 4*x then op([x,y]) else op([0,0]) fi end; r := proc(x, y) if 0 ≤ x and − 2 ≤ y and y ≤ 2 − 4 ∗ x then op([x, y]) else op([0, 0]) end if end proc > f := (x,y)->x*y + x^2; f := (x, y) → x y + x2 @ is the composition operator in Maple. > h :=f @ r;; h := f@r > h(1/3,2/3); 1 3 Here is the graph of f, with the absolute minimum and maximum values shown. > display(pointplot3d({[1,-2,h(1,-2)],[1/3,2/3,h(1/3,2/3)] },color=red,symbol=circle), plot3d(h, > -.1..1.2,-2.1..2.1,view=-1..3/5,style=patchcontour,numpoints=4000,orie ntation=[160,77],axes=boxed)); 10 0 0.5 1 –2–1012 –1 –0.8 –0.6 –0.4 –0.2 0 0.2 0.4 0.6 2. Optimize f(x, y) = x y on the closed disk x2 + y2 ≤ 1 you do this one. Also try to draw a picture to verify your conclusions visually. First change the definition of r so that r(x,y) = (x,y) if xˆ2 + yˆ2 <= 1, otherwise r(x,y) = (0,0) > r:= proc(x,y) if x>= 0 and y >=-2 and y <= 2 - 4*x then op([x,y]) else op([0,0]) fi end; r := proc(x, y) if 0 ≤ x and − 2 ≤ y and y ≤ 2 − 4 ∗ x then op([x, y]) else op([0, 0]) end if end proc the definition of f remains unchanged. > f := (x,y)->x*y + x^2; f := (x, y) → x y + x2 Also, h will remain unchanged as the composition of f with r. > h :=f @ r;; h := f@r Here is the picture you need to modify. You will need to modify the location of the extreme points on the graph. > display(pointplot3d({[1,-2,h(1,-2)],[1/3,2/3,h(1/3,2/3)] },color=red,symbol=circle), plot3d(h, -.1..2,-2.1..2.1,view=-1..1,style=patchcontour,numpoints=4000)); 11 What happens to xbar as T appoaches 0? Since sin(T ) T goes to 1 and 1 − cos(T ) T goes to 0 we see the centroid approaches the point [2R3 , 0]! So, there is a jump discontinuity in the locus at T = 0. Other questions (In all of these, consider R to be fixed.): What is the maximum value of ybar? What is the minimum value of xbar? What is the length of the locus? What is area under the locus? What is the centroid of the region under the locus? Other questions. Take any nice 1 parameter family of regions and ask for the locus of the centroid. For example, it was asked in class where the centroid of a trapezpoidal region is? Making this question very specific, for the right triangle with vertices [0,0], [a,0], [0,b] and h between 0 and b, let T be the trapezoid with base [0,0] to [a,0] and top [0,h] to [a (1 − h b ), h]. Calulate the centroid of T as a function of a, b, and h. Use this to investigate the locus of the centroid with a and b fixed and h going from 0 up to b. Draw pictures to check your work. 14 > restart:with(plottools):with(plots): 8 The centroid of an ice cream cone (spherical box) Fix φ0 between 0 and π. We want to compute the centroid of the solid Sφ0 consisting of all points ( ρ, φ, θ) such that ρ is between 0 and 1, φ is between 0 and φ0, and θ is between 0 and 2π. Here is it’s picture for φ0 = π 4 . > icc:=(phi0,n)->display(plot3d(1, t=0..2*Pi, p=0..phi0, coords=spherical, style=wireframe,color=blue),seq(line([0,0,0],[sin(phi0)*cos(2*Pi*i/n), > sin(phi0)*sin(2*Pi*i/n),cos(phi0)],thickness=2,color=red),i=1..n),scal ing=constrained,axes=normal,orientation=[80,70]): > icc(Pi/4,72); 0.2 0.4 0.6 0.8 1 –0.5 0.5 –0.6–0.4–0.20.20.40.6 > Int(x^3,x=0..3); ∫ 3 0 x3 dx > evalf(%); 20.25000000 By inspection, we can see that the centroid is (0,0,a) where a is somewhere between 0 and 1, depending on φ0. The volume and moment about the xy-plane can best be evaluated as iterated integrals in spherical coordinates. > Int(Int(Int(1,z),x=S[phi[0]]..‘‘),y)=Int(Int(Int(1*rho^2*(sin(phi)),r ho=0..1),phi=0..phi0),theta=0..2*Pi); ∫ ∫ Sφ0 ∫ 1 dz dx dy = ∫ 2 π 0 ∫ φ0 0 ∫ 1 0 ρ2 sin(φ) dρ dφ dθ So the volume of the ice cream cone with cone angle 2φ0 and cone slant height 1 is > V:=int(int(int(1*rho^2*(sin(phi)),rho=0..1),phi=0..phi0),theta=0..2*P i);; V := 2π 3 − 2 3 cos(φ0)π And the moment about the xy plane is 15 > Int(Int(Int(z,z),x=S[phi[0]]..‘‘),y)=Int(Int(Int(rho*cos(phi)*rho^2*si n(phi),rho=0..1), phi=0..phi[0]),theta=0..2*Pi); ∫ ∫ Sφ0 ∫ z dz dx dy = ∫ 2 π 0 ∫ φ0 0 ∫ 1 0 ρ3 cos(φ) sin(φ) dρ dφ dθ This evaluates to > Mxy:=int(int(int(rho*cos(phi)*rho^2*sin(phi),rho=0..1),phi=0..phi0),t heta=0..2*Pi); Mxy := π 4 − 1 4 cos(φ0)2 π So we can calculate the z-coordinate of the centroid, zbar as > zbar:= unapply(Mxy/V,phi0); zbar := φ0 → π 4 − 1 4 cos(φ0)2 π 2π 3 − 2 3 cos(φ0)π And for example when the cone angle is 90 degrees, zbar is > simplify(zbar(Pi/4))=evalf(zbar(Pi/4)); − 3 8 (−2 + √ 2) = 0.6401650420 We can add the centroid to the picture and check to see if it looks reasonable. > picture:= phi0->display(pointplot3d([0,0,evalf(zbar(phi0))],symbol=circle,color= black,thickness=3), icc(phi0,20),scaling=constrained,orientation=[116,80]); picture := φ0 → plots : −display(plots : −pointplot3d ([0, 0, evalf(zbar(φ0))], symbol = plottools : −circle, color = black , thickness = 3), icc(φ0, 20), scaling = constrained , orientation = [116, 80]) > picture(Pi/4); 0.2 0.4 0.6 0.8 1 0.5 –0.6 –0.4–0.2 0.40.6 This seems eminately reasonable. Now lets make a movie of the locus of the centroid of the spherical box as φ0 goes from 0 to Pi. > display(seq(picture(i*Pi/20),i=1..20),scaling=constrained,insequence= true); 16 > frame:= n->display(seq(sq(i*Pi/36,black),i=1..n),inter,inter2,labels=[x,y,z],s tyle=wireframe,orientation=[66,50]); frame := n → plots : −display(seq(sq( i π 36 , black), i = 1..n), inter , inter2 , labels = [x, y, z], style = wireframe , orientation = [66, 50]) > pic:=display(frame(36),axes = boxed): > pic; –1 –0.5 0 0.5 1 x –1 –0.5 0 0.5 1 y –1 –0.5 0 0.5 1 z From this picture we can see that the cross-section by a plane at height h between -1 and 1 is a square of side 2 √ 1 − h2. Hence the volume of the intersection of the two cylinders is ∫ 1 −1 4 − 4h2 dh = 16 3 > Int(4*(1-h^2),h=-1..1) = int(4*(1-h^2),h=-1..1); ∫ 1 −1 4 − 4h2 dh = 16 3 Further questions: What is the surface area of the intersection? also what is the total length of the edges of the intersection. 9.2 Three Cylinders Now lets draw the intersection of the 3 mutually perpendicular cylinders of radius 1 about the x, y, and z axes. The intersection of the cylinder about the z-axis with the other two cylinders in succession give 4 more ellipses which are parameterized and drawn below. > inter3 := display(spacecurve([sin(t),cos(t),sin(t)],t=0..2*Pi,thickness=3,color= black)): > inter4 := display(spacecurve([-sin(t),cos(t),sin(t)],t=0..2*Pi,thickness=3,color =red)): > inter5 := display(spacecurve([cos(t),sin(t),sin(t)],t=0..2*Pi,thickness=3,color= blue)): > inter6 := display(spacecurve([cos(t),-sin(t),sin(t)],t=0..2*Pi,thickness=3,color =brown)): 19 > display(inter5,inter6,inter4,inter3,inter,inter2,labels=[x,y,z],scali ng=constrained,orientation=[40,60],axes=boxed); –1 –0.5 0 0.5 1 x –1 –0.5 0 0.5y –1 –0.5 0 0.5 1 z We can see that the surface of the intersection consists of 12 congruent quadrilateral pieces of a cylinder. They are grouped into 3 groups of 4 pieces of each of the intersecting cylinders. Here is a color coded picture of the surface. > f := proc(t) options operator; if evalf(t)<=evalf(Pi/4) then [sin(t),sin(t),cos(t)] > else [cos(t),sin(t),cos(t)] fi; end: g := proc(t) options operator; > if evalf(t)<=evalf(Pi/4) then [-sin(t),sin(t),cos(t)] else [-cos(t),sin(t),cos(t)] fi; end: > surf:=clr->display(seq(line(f(i*Pi/(2*36)),g(i*Pi/(2*36)),color=clr), i=1..36),seq(pointplot3d({f(i*Pi/(2*36)),g(i*Pi/(2*36)) },symbol=circle,color=black),i=0..36)): > bluesurf:=rotate(display(seq(line(f(i*Pi/(2*36)),g(i*Pi/(2*36)),color =blue),i=1..36)),0,0, Pi/2): > pic:=display(seq(rotate(bluesurf,0,i*Pi/2,0),i=0..3),seq(rotate(refle ct(surf(green),[[0,0,0],[0,0,1],[1,0,0]]),i*Pi/2,0,0 ),i=0..3), seq(rotate(rotate(surf(red),0, Pi/2,0),0,0,i*Pi/2),i=0..3), > scaling=constrained,labels=[x,y,z],axes=boxed,orientation=[70,50],titl e="Intersection of 3 cylinders"): > pic; Intersection of 3 cylinders –1–0.500.51 x –1 –0.5 0 0.5 1 y –1 –0.5 0 0.5 1 z 20 Further examination shows there are 6 vertices of order 4, 2 for each pair of intersecting cylinders and 8 vertices of order 3, representing the 23 = 8 points of intersection of the 3 cylinders. There are also 24 edges. It makes for ball you wouldn’t want to play soccer with. What about the volume? Well, by the symmetry it will 16*V where V is the part trapped above the region 0 <= x <= y and x2 + y2 ≤ 1 in the xy-plane and the cylinder x2 + z2 = 1. This can be written as the sum of two iterated integrals: ∫ 1√ 2 0 ∫ x 0 √ 1 − x2 dy dx + ∫ 1 1√ 2 ∫ √ 1−x2 0 √ 1 − x2 dy dx > V1:=int(sqrt(1-x^2)*x,x=0..1/sqrt(2)); V1 := √ 2 √ π ( 4 √ 2 3 √ π − 2 3 √ π ) 8 > V2:=int(sqrt(1-x^2)*sqrt(1-x^2),x=1/sqrt(2)..1); V2 := 2 3 − 5 √ 2 12 So the total volume of the intersection is > 16*(V1+V2)=evalf(16*(V1+V2)); 2 √ 2 √ π ( 4 √ 2 3 √ π − 2 3 √ π ) + 32 3 − 20 √ 2 3 = 4.686291496 Further questions: Calculate the surface area of the intersection. calculate the total length of the edges of the surface. 21 > onstrained); end use; end: > display(drawit(2,3,Pi/6,4,2,3),translate(drawit(1.5,2,Pi/4,2,2,4),3,0 ,0),orientation=[-40,57],axes=boxed,labels=[x,y,evaln(z)]); 3.535 6.500 0 1 2 3 4 x 00.5 11.5 y 0 1 2 3 4 z Problem: Use this same idea to develop a formula for the volume of a truncated rectangular solid with base rectangle a by b, and three given heights h1, h2, and h3. 24 10 Line integral problems Problem: Let C be the curve sitting on the surface of the graph of smooth function z = f(x,y) over the arc from (0,0) to (1,2), and let g(x,y,z) be a nice smooth density function defined on C. Represent the mass and center of mass of C in terms of line integrals. Evaluate these integrals for f(x,y) = x2 + x y2 and g(x,y,z) = z. Draw a picture showing the surface f and the curve C together with its center of mass. A solution: The mass of C is M = ∫ C g(x, y, z) ds. The moment about the xy-plane is Mx, y = ∫ C z g(x, y, z) ds . The moments about the other two planes are defined in the same way. We can parameterize C by r(t) = (x,y,z) = (t,2*t,f(t,2*t)). So we can evaluate these 4 line integrals like so > M := Int(2*t^2*sqrt(1+4+16*t^2),t=0..1); M := ∫ 1 0 2 t2 √ 5 + 16 t2 dt > M:=int(2*t^2*sqrt(1+4+16*t^2),t=0..1);; M := −5 8 5 √ π 128 − 5 64 ( 1 2 − 6 ln(2) + ln(5))√π − 37 √ π √ 21 40 + 5 32 √ π ln( 1 2 + √ 21 8 ) √ π > M:= evalf(M); M := 2.517952442 > Mxy:=Int((2*t^2)^2*sqrt(1+4+16*t^2),t=0..1); Mxy := ∫ 1 0 4 t4 √ 5 + 16 t2 dt > Mxy:=int((2*t^2)^2*sqrt(1+4+16*t^2),t=0..1); Mxy := −5 4 −125 √ π 12288 + 25 2048 ( 5 6 − 6 ln(2) + ln(5))√π − 711 √ π √ 21 1280 − 25 1024 √ π ln( 1 2 + √ 21 8 ) √ π > Mxy:=evalf(Mxy); Mxy := 3.222893594 > zbar:=Mxy/M; zbar := 1.279966031 > Myz:=Int((2*t^2)*t*sqrt(1+4+16*t^2),t=0..1); Myz := ∫ 1 0 2 t3 √ 5 + 16 t2 dt > Myz:=int((2*t^2)*t*sqrt(1+4+16*t^2),t=0..1); 25 Myz := − 5 (−133 √ π √ 21 200 − √ π √ 5 24 ) 8 √ π > Myz:=evalf(Myz); Myz := 1.962863960 > xbar:=Myz/M; xbar := 0.7795476703 Note that since y = 2*x, Mxz = 2Myz. and so ybar = 2*xbar > Mxz:=2*Myz; Mxz := 3.925727920 > ybar:=Mxz/M; ybar := 1.559095341 Now to draw the picture: > use plots in display(pointplot3d([xbar,ybar,zbar],symbol=circle,color=red), #put this in last > spacecurve([t,2*t,2*t^2],t=0..1,color=blue,thickness=3), #then added this plot3d(x^2+x*y/2,x=0..1,y=0..2,color=yellow), #drew this first > scaling=constrained, #these plot options were added to make the picture more legible. axes=boxed, > labels=[x,y,z], orientation=[147,72]); > end use; 0 0.5 1 x 00.511.5 y 0 0.5 1 1.5 2 h1+h2/a*x-(h1*a+h2*b*cos(t)-h3*a)/b/sin(t)/a*y If this seems a little high, recall that the density function is the z-coordinate and so the center of mass is pulled up significantly toward the upper end of the wire. Problem for you: Recalculate the center of mass if the density of the wire is constant (say 1). Redraw the wire. 26 –1 –0.5 0 0.5 1 y –1 –0.5 0.5 1 x where we note that only the arrows very close to (0,0) are visible. Now using fieldstrength=log > fieldplot([x/(x^2+y^2)^(3/2),y/(x^2+y^2)^(3/2)], x=-1..1,y=-1..1,fieldstrength=log); –1 –0.5 0.5 1 y –1 –0.5 0.5 1 x which makes the direction of the arrows much more visible. Alternatively a pure direction plot can be produced: > fieldplot([x/(x^2+y^2)^(3/2),y/(x^2+y^2)^(3/2)], x=-1..1,y=-1..1,fieldstrength=fixed); –1 –0.5 0 0.5 1 y –1 –0.5 0.5 1 x A radial field in polar coordinates > fieldplot([r,0],r=0..1,t=0..Pi/2,coords=polar); 29 0 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 A field in polar coordinates that will not draw the arrows at the origin (the field direction is undefined there). > fieldplot([0,1],r=0..1,t=0..Pi/2,coords=polar); 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 30 12 Greens theorem. Green’s theorem states that if F = <P,Q> has continuous partial derivatives on an open set containing the positively oriented simple closed curve C and the region D it bounds, then the line integral of the the tangential component of F with repect to arc length around C is the double integral over D of the difference ( ∂ ∂x Q) − ( ∂ ∂y P ). The left hand side of Green’s theorem has a nice meaning: If the integral is positive, then on average the tangential component of F is positive, so if we take the view that F is a velocity flow for a fluid, that says that the net flow of the fluid on the boundary curve C is counterclockwise. If we divide by the length of the curve, we would get average signed speed that a point on the curve is rolled around the curve by the velocity field. If the curve is a small circle of radius r then Green’s theorem says that this average angular velocity is approximately one half the above difference of partial derivatives. For this reason the the difference is a measure of the rotational tendency of the velocity field at each point, and as such is called the curl of F. Note: The curl of a 3 dimensional velocity field F = <P, Q, R> is the vector Del X F = < ( ∂ ∂y R) − ( ∂ ∂z Q), ( ∂ ∂z P ) − ( ∂ ∂x R), ( ∂ ∂y Q) − ( ∂ ∂x P )> . If F is 2 dimensional (R = 0) we can call the 3rd component the curl of F. Stokes theorem is a generalization of Greens theorem to the situation in space where C is a space curve and D is a surface bounded by C. It says that the line integral of the tangential component of F around C is equal to the surface integral of the normal component of curl F over D. If instead of integrating the tangential component of the velocity field F around the curve C, we integrate the normal component of F, we get another line integral which can also be equated to a double integral by Green’s theorem. The outward normal n to the curve C is obtained by rotating the tangent vector <dx/dt, dy/dt>/—r’(t)— 90 degrees clockwise to get n = <dy/dt, -dx/dt>/—r’(t)— The line integral of the normal component of F = <P,Q> then becomes ∫ C P dy − ∫ C Qdx. By Greens theorem this is equal to the the double integral ∫ ∫ ( ∂ ∂x P ) + ( ∂ ∂y Q) d dA, This is called the normal form of Greens theorem. The right hand side of this equation measures the transport of fliud across the curve C (positive means fluid is flowing out on average, negative means fluid is flowing in). The integrand of the right hand integral is called the divergence of the velocity field F. More generally, the divergence of a vector field F = <P,Q,R> is defined as div F = del dot F = ( ∂ ∂x P ) + ( ∂ ∂y Q) + ( ∂ ∂z R). The Divergence theorem is a generalization of the normal form of Green’s theorem to the situation in space where you have a surface R bounding a soliid S in space where there is a nice velocity field F defined on an open set containing S and R. It says that the surface integral of the (outward) normal component of F over R is equal to the triple integral of the divergence of F over the solid S. 12.1 procedure 1 to investigate velocity fields: Tangential form of Greens the- orem. This procedure takes a vector field F1 (a list of two expressions in x and y), a radius R, and a center P and draws the circle C of radius R centered at P. Then it draws the the red graph of the 31