1.14 Lesson 14

```# 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
> 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);

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

```