# Lesson 14. (a) Optimization on a Restricted Domain (b) Introduction to Integration # ----------------------------------------------------------------------------- # ------- # # Suppose we wish to find the maximum of a function f(x,y) on a domain D. The maximum (if it exists) # could occur at an interior point of D, in which case it is a critical point or singular point of f. Or it could # occur at a boundary point of D. Typically the boundary of D might consist of one or more smooth # curves, each defined by an equation g[j](x,y) = 0. A maximum on one of these curves might be found # using Lagrange multipliers. Or the maximum could occur at the intersection of two boundary curves. # All these possibilities must be investigated. For a function of three variables the situation is even more # complicated: the boundary consists of surfaces, which intersect on lines, which in turn intersect at points. > restart; with(plots): with(linalg): Warning: new definition for norm Warning: new definition for trace # Example: Find the maximum and minimum of x^2 - xy + y^3 - x # subject to the constraints x^2 + x y + y^2 <= 3 and x^2 + y^2 <= 4. > f:= (x,y) -> x^2 -x*y + y^3 - x: > g1:= (x,y) -> x^2 + x*y + y^2 - 3: > g2:= (x,y) -> x^2 + y^2 - 4: # Let's first try a plot, to see what the domain and the objective function x^2 - xy + y^3 - x look like, and # make an educated guess at the locations of the maximum and minimum. What type of plots shall we use? # An "implicitplot" would be a good choice to show the boundary of our domain. A "contourplot", on the # other hand, might be good to show the objective. In Release 4 we could combine these with "display". In # Release 3, unfortunately, "implicitplot" produces a two-dimensional plot while "contourplot" produces a # three-dimensional one, and these can't be combined. Instead, we can simulate a contour plot using # "implicitplot". > plot1:= implicitplot({g1(x,y)=0,g2(x,y)=0}, > x=-2..2, y=-2..2, thickness=3): > plot2:= implicitplot({seq(f(x,y)=3*kk, kk=-2..3)}, > x=-2..2, y=-2..2): > display({plot1,plot2}, axes=BOXED); # Our domain is the region inside both thick curves. The contour values increase # from bottom to top (note, for example, that f(0,y) = y^3 is an increasing function # of y). According to the plot, it appears that the maximum value, approximately # 9, occurs on the curve g2(x,y) = 0 somewhere between about (-1,1.7) and the intersection of the two # curves near (-.5, 1.95) (the curve g2(x,y)=0 and the contour f(x,y) = 9 are almost indistinguishable in this # region), while the minimum # value, slightly less than -6, occurs at the intersection of the two curves near (.5, -1.95). # # Now we'll use calculus. There are four places to search: # 1) Critical points of f in the region g1(x,y) < 0, g2(x,y) < 0. # 2) Critical points of the restriction of f to g1(x,y) = 0 with g2(x,y) < 0. # 3) Critical points of the restriction of f to g2(x,y) = 0 with g1(x,y) < 0. # 4) Intersections of g1(x,y)=0 and g2(x,y) = 0. # # I'll just do (2) here, and you can do the rest in your homework. # # We introduce Lagrange multipliers. > L:= (x,y,t) -> f(x,y) + t*g1(x,y): > eqs:= { diff(L(x,y,t),x)=0, diff(L(x,y,t),y)=0, > diff(L(x,y,t),t)=0 }; 2 2 eqs := {x + x y + y - 3 = 0, 2 x - y - 1 + t (2 x + y) = 0, 2 - x + 3 y + t (x + 2 y) = 0} > sols:= solve(eqs, {x,y,t}); 72 4 25 3 34 2 901 2719 321 5 sols := {t = - ---- %1 + --- %1 - ---- %1 - ---- + ---- %1 - ---- %1 , 1355 271 1355 1355 4065 1355 54 5 468 4 27 3 1134 2 658 1596 y = %1, x = - ---- %1 - ---- %1 + --- %1 + ---- %1 - ---- %1 + ----} 1355 1355 271 1355 1355 1355 6 4 3 2 %1 := RootOf(9 _Z - 21 _Z + 6 _Z - 59 _Z - 12 _Z + 47) # The equations being polynomials, these should be all solutions. > ro:= subs(sols, y); 6 4 3 2 ro := RootOf(9 _Z - 21 _Z + 6 _Z - 59 _Z - 12 _Z + 47) > ys:= [fsolve(op(ro), _Z)]; ys := [-1.961526843, -.8752838891, .7701014499, 1.896406235] > cps:= map(r -> subs(ro=r, sols), ys); cps := [{t = 3.322808021, y = -1.961526843, x = .6426667844}, {t = -1.240972490, y = -.8752838891, x = 1.995014278}, {x = 1.213451017, y = .7701014499, t = -.2054425610}, {t = -5.355428469, y = 1.896406235, x = -1.498414482}] # Check for accuracy. > map(S -> subs(S, eqs), cps); -8 -7 -7 -9 [{.3*10 = 0, -.2*10 = 0, -.23*10 = 0}, {0 = 0, .2*10 = 0}, -9 -8 -9 -8 -7 {.7*10 = 0, .2*10 = 0, -.8*10 = 0}, {0 = 0, .2*10 = 0, .12*10 = 0} ] # Check to see if these are in our region (g2(x,y) < 0). > map(S -> subs(S, g2(x,y)), cps); [.260608152, .746203856, -1.934480386, 1.841602568] > val:= subs(cps[3], f(x,y)); val := -.2187545609 # # Integration # ----------- # # The main integration command in Maple is "int". It can be used for indefinite integrals: > int(x^2, x); 3 1/3 x # or definite integrals: > int(x^2, x = 1 .. 4); 21 # Double and triple integrals can be done as iterated integrals: > int(int(x^2*y, x = 1 .. 2), y = 3 .. 4); 49/6 # There is also an "inert" integration command Int for when you don't really want to do the integration. > Int(x^2, x = 1 .. 4); 4 / | 2 | x dx | / 1 # The result of Int can be manipulated in various ways. When you want to evaluate the integral, you can # use "value" on the expression containing "Int". > " = value("); 4 / | 2 | x dx = 21 | / 1 # As you know, it can be hard to find antiderivatives. Maple is pretty good at it, but not perfect. Here's a # hard one it can do. > J:=int(sqrt(tan(x)), x); 1/2 1/2 1/2 1/2 1/2 2 tan(x) 1/2 tan(x) + 2 tan(x) + 1 J := 1/2 2 arctan(--------------) - 1/2 2 ln(---------------------------) 1 - tan(x) 2 1/2 (1 + tan(x) ) # The way to check an antiderivative is to differentiate it. > dJ:= diff(J, x); / 1/2 2 1/2 1/2 2 \ 1/2 | 2 (1 + tan(x) ) 2 tan(x) (- 1 - tan(x) )| 2 |1/2 ---------------------- - ------------------------------| | 1/2 2 | \ tan(x) (1 - tan(x)) (1 - tan(x)) / dJ := 1/2 ------------------------------------------------------------------ tan(x) 1 + 2 ------------- 2 (1 - tan(x)) 1/2 - 1/2 2 / 1/2 2 \ | 2 2 (1 + tan(x) ) | |1 + tan(x) + 1/2 ------------------ | | 1/2 1/2 1/2 | | tan(x) (tan(x) + 2 tan(x) + 1) tan(x)| |------------------------------------ - ------------------------------------| | 2 1/2 2 1/2 | \ (1 + tan(x) ) (1 + tan(x) ) / 2 1/2 / 1/2 1/2 (1 + tan(x) ) / (tan(x) + 2 tan(x) + 1) / # Doesn't look much like the original. Maybe it simplifies. "normal" does a moderate level of # simplification (putting everything over a common denominator, taking out factors). # "simplify" would also do other things, such as replacing tan(x) by sin(x)/cos(x), which we don't want. > normal("); 1/2 1/2 2 3/2 1/2 2 (2 tan(x) + 2 tan(x) + 2 tan(x)) 1/2 ----------------------------------------------- 1/2 1/2 1/2 tan(x) (tan(x) + 2 tan(x) + 1) # You can probably recognize this as equivalent to sqrt(tan(x)) if you stare at it long enough. It's hard to # get Maple to simplify it to sqrt(tan(x)). To check whether a = b it's often better to try checking whether a # - b = 0. > normal(dJ - sqrt(tan(x))); 0 # A harder one: > int(sqrt(sec(x)^2 + 4), x); # We can give Maple a little help here: try a change of variables. # # It turns out that the substitution t = tan(x) works. You can't # simply use "subs" for a change of variables in an integral because # you have to substitute in the "dx" as well. The "student" package # has a command called "changevar" that will do it for us. # # Up to now, when we wanted to use a command from a package, # we used "with". You don't really need to do that: you can use # the name of the package followed by the name of the command # in brackets. If you're going to use the command more than once, # or use a number of different commands from the package, "with" # is still worthwhile. > student[changevar](tan(x)=t, ", t); t 1/2 2 arctan(2 -----------) + arcsinh(1/5 5 t) 2 1/2 (5 + t ) > subs(t=tan(x),"); tan(x) 1/2 2 arctan(2 ----------------) + arcsinh(1/5 5 tan(x)) 2 1/2 (5 + tan(x) ) > normal(diff(",x) - sqrt(sec(x)^2+4)); 2 1/2 2 1/2 2 - (- 20 (25 + 5 tan(x) ) - 4 (25 + 5 tan(x) ) tan(x) 1/2 2 3/2 1/2 2 3/2 2 - 5 (5 + tan(x) ) - 5 (5 + tan(x) ) tan(x) 2 1/2 2 3/2 2 1/2 + (sec(x) + 4) (5 + tan(x) ) (25 + 5 tan(x) ) ) / 2 3/2 2 1/2 / ((5 + tan(x) ) (25 + 5 tan(x) ) ) / > simplify("); 0 # Some antiderivatives can't be written as "elementary functions" (i.e. combinations of the functions Math # 100 students know about). > int(exp(sin(x)), x); / | | exp(sin(x)) dx | / # Some of these antiderivatives can be expressed in terms of special functions that Maple knows about. # Here's one that's important in probability. > int(exp(-x^2), x); 1/2 1/2 Pi erf(x) # Here's another that you may not know about. > int(ln(x)/(1+x), x); dilog(1 + x) + ln(x) ln(1 + x) > ?dilog # FUNCTION: dilog - The Dilogarithm Function # # CALLING SEQUENCE: # dilog(x) # # PARAMETERS: # x - an expression # # SYNOPSIS: # - The dilogarithm function is defined as follows: # # dilog(x) = int(ln(t)/(1-t), t=1..x) > int(ln(x)/(1+x), x= 0 .. 1); 2 - 1/12 Pi # When Maple can find an antiderivative of a function, it should have no trouble with definite integrals of # that function, using the Fundamental Theorem of Calculus. Well, maybe... > int(sqrt(tan(x)), x = 0 .. Pi/3); 1/2 1/4 1/2 2 3 1/2 1/2 1/2 1/4 1/2 - 1/2 2 arctan(----------) - 1/2 2 ln(1/4 (3 + 2 3 + 1) 4 ) 1/2 - 1 + 3 1/2 + 1/2 2 Pi > evalf("); .787779048 # Knowing the antiderivative (which we called J), could we have # done that ourselves? > subs(x=Pi/3, J) - subs(x=0, J); 1/2 1/2 1/2 2 tan(1/3 Pi) 1/2 2 arctan(-------------------) 1 - tan(1/3 Pi) 1/2 1/2 1/2 tan(1/3 Pi) + 2 tan(1/3 Pi) + 1 - 1/2 2 ln(-------------------------------------) 2 1/2 (1 + tan(1/3 Pi) ) 1/2 1/2 1/2 2 tan(0) - 1/2 2 arctan(--------------) 1 - tan(0) 1/2 1/2 1/2 tan(0) + 2 tan(0) + 1 + 1/2 2 ln(---------------------------) 2 1/2 (1 + tan(0) ) > evalf("); -1.433662421 # Oops! What went wrong? # Which answer is more reasonable? # Take another look at J. > J; 1/2 1/2 1/2 1/2 1/2 2 tan(x) 1/2 tan(x) + 2 tan(x) + 1 1/2 2 arctan(--------------) - 1/2 2 ln(---------------------------) 1 - tan(x) 2 1/2 (1 + tan(x) ) > plot(J, x=0 .. Pi/3); > readlib(discont): > discont(J, x); {Pi _Z7~ + 1/2 Pi, 1/4 Pi + Pi _Z~} # What about if we had an example where Maple couldn't identify # the discontinuities? > f:= 1/(2*x - tan(x)); 1 f := ------------ 2 x - tan(x) > g:= normal(diff(f,x)); 2 - 1 + tan(x) g := --------------- 2 (2 x - tan(x)) > J:=int(g,x); 1 - -------------- - 2 x + tan(x) > evalf([subs(x=Pi/4, J), subs(x=3*Pi/8, J)]); [1.751938393, -17.23571203] > int(g, x=Pi/4 .. 3*Pi/8); 4 2 ----------------- - ------ 1/2 Pi - 2 3 Pi - 4 2 - 4 > evalf("); -18.98765027 # What should the answer be?