# 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?
