# Solutions to Homework 4 (Sjerve)

```#
#
#                                 HOMEWORK ASSIGNMENT # 4,
#
#                                           Due March 15, 1996
#
#
--------------------------------------------------------------------------------
#  QUESTION #1
#
#      Find all solutions of the following system of equations
#
#                         x+y+z=-2, x^2+y^2+z^2=-2, x^3+y^3+z^3=-2.
#
--------------------------------------------------------------------------------
# SOLUTION:
--------------------------------------------------------------------------------
> restart;
--------------------------------------------------------------------------------
> eq1:=x+y+z=-2;eq2:=x^2+y^2+z^2=-2;eq3:=x^3+y^3+z^3=-2;

eq1 := x + y + z = -2

2    2    2
eq2 := x  + y  + z  = -2

3    3    3
eq3 := x  + y  + z  = -2
--------------------------------------------------------------------------------
# From eq1 we see that z=-x-y-2, so we substitute this into eq2 and eq3, obtaining 2
# equations in x and y.
--------------------------------------------------------------------------------
> eq4:=subs(z=-x-y-2,eq2);

2    2                2
eq4 := x  + y  + (- x - y - 2)  = -2
--------------------------------------------------------------------------------
> eq5:=subs(z=-x-y-2,eq3);

3    3                3
eq5 := x  + y  + (- x - y - 2)  = -2
--------------------------------------------------------------------------------
> f:=x^2+y^2+(x+y+2)^2;

2    2              2
f := x  + y  + (x + y + 2)
--------------------------------------------------------------------------------
> g:=x^3+y^3-(x+y+2)^3;\

3    3              3
g := x  + y  - (x + y + 2)
--------------------------------------------------------------------------------
# We want to find all (x,y) such that f(x,y)=-2 and g(x,y)=-2.  Then the solutions to the
# original equations are (x,y,z)=(x,y,-x-y-2).  Thus we first solve for y in the quadratic
# equation f(x,y)=-2, obtaining 2 functions of x, which we then substitute into the cubic
# equation g(x,y)=-2.  We expect a total of 6 solutions.
--------------------------------------------------------------------------------
> Y:=solve(f=-2,y);

Y :=

2           1/2                          2           1/2
- 1/2 x - 1 + 1/2 (- 3 x  - 4 x - 8)   , - 1/2 x - 1 - 1/2 (- 3 x  - 4 x - 8)
--------------------------------------------------------------------------------
> Y1:=unapply(Y[1],x);

2           1/2
Y1 := x -> - 1/2 x - 1 + 1/2 (- 3 x  - 4 x - 8)
--------------------------------------------------------------------------------
> Y2:=unapply(Y[2],x);

2           1/2
Y2 := x -> - 1/2 x - 1 - 1/2 (- 3 x  - 4 x - 8)
--------------------------------------------------------------------------------
> h(x):=expand(subs(y=Y1(x),g));

3            2
h(x) := 3 x  + 9 x + 6 x  + 10
--------------------------------------------------------------------------------
# Substituting y=Y2(x) into g will produce the same cubic, so we only need h(x).
--------------------------------------------------------------------------------
> h1(x):=expand(subs(y=Y2(x),g));

3            2
h1(x) := 3 x  + 9 x + 6 x  + 10
--------------------------------------------------------------------------------
> X:=solve(h(x)=-2,x);

5
X := - %2 + ---------------------- - 2/3,
/ 35         1/2\1/3
9 |---- + 5/9 6   |
\ 27            /

5
1/2 %2 - ----------------------- - 2/3
/ 35         1/2\1/3
18 |---- + 5/9 6   |
\ 27            /

1/2 /                  5          \
+ 1/2 I 3    |- %2 - ----------------------|,
|         / 35         1/2\1/3|
|       9 |---- + 5/9 6   |   |
\         \ 27            /   /

5
1/2 %2 - ----------------------- - 2/3
/ 35         1/2\1/3
18 |---- + 5/9 6   |
\ 27            /

1/2 /                  5          \
- 1/2 I 3    |- %2 - ----------------------|
|         / 35         1/2\1/3|
|       9 |---- + 5/9 6   |   |
\         \ 27            /   /

1
%1 :=                        --------------------
/ 35         1/2\1/3
|---- + 5/9 6   |
\ 27            /

/ 35         1/2\1/3
%2 :=                        |---- + 5/9 6   |
\ 27            /
--------------------------------------------------------------------------------
# These expressions are the x coordinates at the points of intersection.
--------------------------------------------------------------------------------
> evalf(X);

-1.650629192,  - .1746854042 - 1.546868888 I,  - .1746854042 + 1.546868888 I
--------------------------------------------------------------------------------
>  for j from 1 to 3 do x[j]:=evalf(X[j]):y[j]:=Y1(evalf(X[j])) od ;j:='j':

x[1] := -1.650629192

y[1] :=  - .1746854040 + 1.546868888 I

x[2] :=  - .1746854042 - 1.546868888 I

y[2] :=  - .1746854034 + 1.546868888 I

x[3] :=  - .1746854042 + 1.546868888 I

y[3] :=  - .1746854034 - 1.546868888 I
--------------------------------------------------------------------------------
>  for j from 1 to 3 do x[j+3]:=evalf(X[j]);y[j+3]:=Y2(evalf(X[j])) od ;j:='j':

x[4] := -1.650629192

y[4] :=  - .1746854040 - 1.546868888 I

x[5] :=  - .1746854042 - 1.546868888 I

-9
y[5] :=  - 1.650629192 + .5*10   I

x[6] :=  - .1746854042 + 1.546868888 I

-9
y[6] :=  - 1.650629192 - .5*10   I
--------------------------------------------------------------------------------
# Now we print out the 6 solutions.
--------------------------------------------------------------------------------
> for j from 1 to 6 do print(S[j]=[x[j],y[j],-x[j]-y[j]-2]) od;j:='j':

S[1] =

[-1.650629192,  - .1746854040 + 1.546868888 I,  - .174685404 - 1.546868888 I]

S[2] =

[ - .1746854042 - 1.546868888 I,  - .1746854034 + 1.546868888 I, -1.650629192]

S[3] =

[ - .1746854042 + 1.546868888 I,  - .1746854034 - 1.546868888 I, -1.650629192]

S[4] =

[-1.650629192,  - .1746854040 - 1.546868888 I,  - .174685404 + 1.546868888 I]

-9
S[5] = [ - .1746854042 - 1.546868888 I,  - 1.650629192 + .5*10   I,

- .174685404 + 1.546868888 I]

-9
S[6] = [ - .1746854042 + 1.546868888 I,  - 1.650629192 - .5*10   I,

- .174685404 - 1.546868888 I]
--------------------------------------------------------------------------------
# Next we check the residuals to make sure all values are close to 2.
--------------------------------------------------------------------------------
> for j from 1 to 6 do subs(x=x[j],y=y[j],f);subs(x=x[j],y=y[j],g) od;j:='j':

-2.000000003

-2.000000006

-8
- 2.000000000 + .25*10   I

-8
- 2.000000018 - .1*10   I

-8
- 2.000000000 - .25*10   I

-8
- 2.000000018 + .1*10   I

-2.000000003

-2.000000006

-8
- 2.000000003 - .1050629192*10   I

-8
- 2.000000005 + .4086865094*10   I

-8
- 2.000000003 + .1050629192*10   I

-8
- 2.000000005 - .4086865094*10   I
--------------------------------------------------------------------------------
> f:='f':g:='g':Y:='Y':Y1:='Y1':Y2:='Y2':h:='h':X:='X':S:='S':
--------------------------------------------------------------------------------
# There is another approach to this question.  Notice that the original equations are
# symmetric in x,y and z, so if we find one solution then we get 5 more by permuting the
# coordinates.
--------------------------------------------------------------------------------
# QUESTION #2
--------------------------------------------------------------------------------
# SOLUTION
--------------------------------------------------------------------------------
# 2(a)
--------------------------------------------------------------------------------
> f:=x->(x+1)/(x^3+x+1);

x + 1
f := x -> ----------
3
x  + x + 1
--------------------------------------------------------------------------------
> F(x):=int(f(x),x);

-----
\                93    2   124       14
F(x) :=    )    _R ln(x + ---- _R  + --- _R + ----)
/                47         47       47
-----
_R = %1

3
%1 :=                     RootOf(31 _Z  + 7 _Z - 1)
--------------------------------------------------------------------------------
# 2(b)
--------------------------------------------------------------------------------
> diff(F(x),x);

x + 1
----------
3
x  + x + 1
--------------------------------------------------------------------------------
> f:='f':F:='F':
--------------------------------------------------------------------------------
# QUESTION #3
--------------------------------------------------------------------------------
# 3(a)
--------------------------------------------------------------------------------
> g:=x->floor(x/Pi);

x
g := x -> floor(----)
Pi
--------------------------------------------------------------------------------
# This function gives the correct values on the interval 0<x<2*Pi.
--------------------------------------------------------------------------------
> plot(g(x),x=0..4*Pi,title=`the graph of floor(x/2*Pi)`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
# The vertical line segments are not part of the graph.  It consists of horizontal line
# segments of length Pi starting at x=0 and jumping up by 1 unit at every multiple of Pi.
--------------------------------------------------------------------------------
#
> gper:=x->g(x-2*Pi*floor(x/(2*Pi)));

x
gper := x -> g(x - 2 Pi floor(1/2 ----))
Pi
--------------------------------------------------------------------------------
> with(plots):
--------------------------------------------------------------------------------
> G:=plot(gper(x),x=-2*Pi..2*Pi,y=-1..2):
--------------------------------------------------------------------------------
> display(G,title=`periodic extension of g(x)`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

#
--------------------------------------------------------------------------------
# 3(b)
--------------------------------------------------------------------------------
> A:=(1/Pi)*int(g(x)*cos(n*x),x=0..2*Pi);

2 Pi
/
|            x
|    floor(----) cos(n x) dx
|           Pi
/
0
A := -----------------------------
Pi
--------------------------------------------------------------------------------
# Since g(x)=0 for 0<x<Pi and g(x)=1 for Pi<x<2*Pi we have, at least for n>0:
--------------------------------------------------------------------------------
> A:=(1/Pi)*int(cos(n*x),x=Pi..2*Pi);

sin(2 n Pi)   sin(n Pi)
----------- - ---------
n            n
A := -----------------------
Pi
--------------------------------------------------------------------------------
> weknow1:=sin(2*n*Pi)=0,sin(n*Pi)=0;

weknow1 := sin(2 n Pi) = 0, sin(n Pi) = 0
--------------------------------------------------------------------------------
> subs(weknow1,A);

0
--------------------------------------------------------------------------------
# Thus all a[n]=0, n>0.  For a[0] we have:
--------------------------------------------------------------------------------
> a[0]=(1/Pi)*int(1,x=Pi..2*Pi);

a[0] = 1
--------------------------------------------------------------------------------
# Notice that gper(x)-1/2 is an odd function.
--------------------------------------------------------------------------------
> B:=(1/Pi)*int(sin(n*x),x=Pi..2*Pi);

cos(2 n Pi)   cos(n Pi)
- ----------- + ---------
n            n
B := -------------------------
Pi
--------------------------------------------------------------------------------
> weknow2:=cos(2*n*Pi)=1,cos(n*Pi)=(-1)^n;

n
weknow2 := cos(2 n Pi) = 1, cos(n Pi) = (-1)
--------------------------------------------------------------------------------
> subs(weknow2,B);

n
(-1)
- 1/n + -----
n
-------------
Pi
--------------------------------------------------------------------------------
# Thus the b[n] coefficients are:
--------------------------------------------------------------------------------
> b[2*n]=0;b[2*n-1]=-2/((2*n-1)*Pi);

b[2 n] = 0

2
b[2 n - 1] = - ------------
(2 n - 1) Pi
--------------------------------------------------------------------------------
# The Fourier series is given by:
--------------------------------------------------------------------------------
> Fourier(x)=1/2-(2/Pi)*sum((1/(2*n-1))*sin((2*n-1)*x),n=1..infinity);

infinity
-----
\      sin((2 n - 1) x)
)     ----------------
/           2 n - 1
-----
n = 1
Fourier(x) = 1/2 - 2 -------------------------
Pi
--------------------------------------------------------------------------------
# 3(c)
--------------------------------------------------------------------------------
> S:=n->1/2-(2/Pi)*sum((1/(2*k-1))*sin((2*k-1)*x),k=1..n);

n
-----
\    sin((2 k - 1) x)
)   ----------------
/         2 k - 1
-----
k = 1
S := n -> 1/2 - 2 ----------------------
Pi
--------------------------------------------------------------------------------
> S(1);S(2);S(3);

sin(x)
1/2 - 2 ------
Pi

sin(x) + 1/3 sin(3 x)
1/2 - 2 ---------------------
Pi

sin(x) + 1/3 sin(3 x) + 1/5 sin(5 x)
1/2 - 2 ------------------------------------
Pi
--------------------------------------------------------------------------------
> for n from 1 to 30 do P.n:=plot(S(n),x=0..2*Pi) od:n:='n':
--------------------------------------------------------------------------------
> Gr:=plot(gper(x),x=0..2*Pi):
--------------------------------------------------------------------------------
> display([Gr,P1],title=`first partial sum`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

#
--------------------------------------------------------------------------------
> display([Gr,P1,P2],title=`first 2 partial sums`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
> display([Gr,P1,P2,P3],title=`first 3 partial sums`);
#
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
> display([Gr,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10],title=`first 10 partial sums`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
> display([Gr,P10,P20,P30],title=`more partial sums`);
--------------------------------------------------------------------------------
#

** Maple V Graphics **

--------------------------------------------------------------------------------
# 3(d)
--------------------------------------------------------------------------------
# For Gibbs' phenomenum we want to find the smallest positive critical point x[n] of S(n)
# and then compute limit(S[n](x[n]),infinity).  Thus we are approximating the excess by
# approaching 0 from the right.  We get similar results if we approach from the left or if we
# approach Pi.
--------------------------------------------------------------------------------
> S:=n->1/2-(2/Pi)*sum((1/(2*k-1))*sin((2*k-1)*x),k=1..n);

n
-----
\    sin((2 k - 1) x)
)   ----------------
/         2 k - 1
-----
k = 1
S := n -> 1/2 - 2 ----------------------
Pi
--------------------------------------------------------------------------------
# The derivative of S(n)(x) is given by f(x) as follows.
--------------------------------------------------------------------------------
> f:=n->(-2/Pi)*sum(cos((2*k-1)*x),k=1..n);

n
-----
\
)   cos((2 k - 1) x)
/
-----
k = 1
f := n -> - 2 ----------------------
Pi
--------------------------------------------------------------------------------
> f(n);

- 2 (- cos(x (n + 1))

2
(2 cos(x) cos(x (n + 1)) sin(x) - 2 cos(x)  sin(x (n + 1)) + sin(x (n + 1)))

/sin(x) + cos(x))/Pi
--------------------------------------------------------------------------------
# Maple knows some trig.
--------------------------------------------------------------------------------
> expand(");

2                             4
sin(x) cos(n x) cos(x)  sin(n x)     cos(n x) cos(x)  sin(n x)
- 4 -------------------------------- - 4 -------------------------
Pi                            Pi sin(x)

2                      2
cos(n x) cos(x)  sin(n x)     cos(n x)  cos(x)
+ 2 ------------------------- + 2 ----------------
Pi sin(x)                    Pi

2         2                    2       3             2
sin(x)  sin(n x)  cos(x)     sin(n x)  cos(x)      sin(n x)  cos(x)
+ 4 ------------------------ + 4 ----------------- - 2 ----------------
Pi                        Pi                   Pi

sin(x) sin(n x) cos(n x)     cos(x)
- 2 ------------------------ - 2 ------
Pi                  Pi
--------------------------------------------------------------------------------
# We want to find positive, but close to zero, values of x so that this last expression is 0.
# For values of x close to 0 we may replace sin(n*x)/sin(x) by n and cos(x) by 1.  We may
# also omit any term that has sin(x) in the numerator and is not zero in the denominator.
# If we do this then we get the following function.
--------------------------------------------------------------------------------
> h:=x->(-4*n/Pi)*cos(n*x)+(2*n/Pi)*cos(n*x)+(2/Pi)*cos(n*x)^2+(2/Pi)*(sin(n*x))^2-2/Pi;

2             2
n cos(n x)     cos(n x)      sin(n x)      2
h := x -> - 2 ---------- + 2 --------- + 2 --------- - ----
Pi             Pi            Pi       Pi
--------------------------------------------------------------------------------
# This simplifies to:
--------------------------------------------------------------------------------
> h:=x->-(2*n/Pi)*cos(n*x);

n cos(n x)
h := x -> - 2 ----------
Pi
--------------------------------------------------------------------------------
> h(0);

n
- 2 ----
Pi
--------------------------------------------------------------------------------
> h(Pi/(2*n));

0
--------------------------------------------------------------------------------
> h(Pi/n);

n
2 ----
Pi
--------------------------------------------------------------------------------
> x[n]:=Pi/(2*n);

Pi
x[n] := 1/2 ----
n
--------------------------------------------------------------------------------
# Thus as x went from 0 to Pi/n, h(x) went from the large negative value -2*n/Pi to the
# large positive value 2*n/Pi, with h(x)=0 at the mid point Pi/(2*n).  Thus we should expect
# f(n) to have a zero quite close to x[n].
--------------------------------------------------------------------------------
> G:=n->evalf(subs(x=Pi/(2*n),S(n))):
--------------------------------------------------------------------------------
> starttime:=time():
--------------------------------------------------------------------------------
> for n from 1 to 500 do G(n) od:n:='n':
--------------------------------------------------------------------------------
> endtime:=time()-starttime;

endtime := 376.650
--------------------------------------------------------------------------------
> G(50);G(60);G(70);G(80);G(90);G(100);

-.0895065396

-.0895014468

-.0894983760

-.0894963828

-.0894950162

-.0894940390
--------------------------------------------------------------------------------
> for n from 480 to 500 do G(n) od;n:='n':

-.0894900540

-.0894900522

-.0894900516

-.0894900510

-.0894900504

-.0894900494

-.0894900490

-.0894900480

-.0894900468

-.0894900464

-.0894900456

-.0894900442

-.0894900442

-.0894900434

-.0894900428

-.0894900422

-.0894900420

-.0894900412

-.0894900398

-.0894900390

-.0894900384
--------------------------------------------------------------------------------
# It would seem that limit(S(n)(x[n],n=infinity) is approximately.-0.08949
--------------------------------------------------------------------------------
# QUESTION #4
--------------------------------------------------------------------------------
>  f:=(x,y)->x^3+y^3-3*x*y;\

3    3
f := (x,y) -> x  + y  - 3 x y
--------------------------------------------------------------------------------
>  f1:=D[1](f);f2:=D[2](f);

2
f1 := (x,y) -> 3 x  - 3 y

2
f2 := (x,y) -> 3 y  - 3 x
--------------------------------------------------------------------------------
> solve({f1(x,y),f2(x,y)},{x,y});\

{y = 0, x = 0}, {y = 1, x = 1},

2                           2
{y = RootOf(_Z  + _Z + 1), x = - RootOf(_Z  + _Z + 1) - 1}
--------------------------------------------------------------------------------
# Thus there are 2 critical points, (x,y)=(0,0), (1,1).  The other solution  is complex.
--------------------------------------------------------------------------------
> with(linalg):\

Warning: new definition for   norm
Warning: new definition for   trace

--------------------------------------------------------------------------------
>  H:=hessian(f(x,y),[x,y]);

[ 6 x   -3 ]
H := [          ]
[  -3  6 y ]
--------------------------------------------------------------------------------
> x:=0;y:=0;

x := 0

y := 0
--------------------------------------------------------------------------------
>  A1:=matrix(2,2,[H[1,1],H[1,2],H[2,1],H[2,2]]);

[  0  -3 ]
A1 := [        ]
[ -3   0 ]
--------------------------------------------------------------------------------
> eigenvals(A1);

3, -3
--------------------------------------------------------------------------------
# The eigenvalues have opposite sign and therefore (0,0) is a saddle point.
--------------------------------------------------------------------------------
>  x:='x';y:='y';

x := x

y := y
--------------------------------------------------------------------------------
>  x:=1;y:=1;

x := 1

y := 1
--------------------------------------------------------------------------------
> A2:=matrix(2,2,[H[1,1],H[1,2],H[2,1],H[2,2]]);

[  6  -3 ]
A2 := [        ]
[ -3   6 ]
--------------------------------------------------------------------------------
>  eigenvals(A2);

9, 3
--------------------------------------------------------------------------------
# The eigenvalues are both positive and therefore (1,1) is a local minimum.
--------------------------------------------------------------------------------
# QUESTION #5
--------------------------------------------------------------------------------
# 5(a)
--------------------------------------------------------------------------------
> int(1/(sin(x)-1/2),x=0..Pi);

1/2           1/2              1/2
- 2/3 I 3    Pi + 4/3 3    arctanh(2/3 3   )
--------------------------------------------------------------------------------
> evalf(");

3.041383981 - 7.255197458 I
--------------------------------------------------------------------------------
# Maple has an answer for this integral, but clearly it is not correct.  The answer is
# complex, but the integrand is real and therefore so should the answer be real.  The trouble
# is that the integrand f(x)=1/(sin(x)-1/2) has singularities at x=Pi/6 and 5*Pi/6, and Maple
# has outsmarted itself in the way it handled these singularities.  The integral actually
# diverges.
--------------------------------------------------------------------------------
> sin(Pi/6);sin(5*Pi/6);

1/2

1/2
--------------------------------------------------------------------------------
> limit((sin(x)-1/2)/(x-Pi/6),x=Pi/6);limit((sin(x)-1/2)/(x-5*Pi/6),x=5*Pi/6);

1/2
1/2 3

1/2
- 1/2 3
--------------------------------------------------------------------------------
# Since the integrals int(1/(x-Pi/6,x=0..Pi) and  int(1/(x-5*Pi/6,x=0..Pi) diverge so does
# int(1/(sin(x)-1/2),x=0..Pi).
--------------------------------------------------------------------------------
# 5(b)
--------------------------------------------------------------------------------
# The integral int(1/(sin(x)-1/2),x=Pi/4..Pi) still diverges because of the singularity at