#
#
# 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
# x=5*Pi/6. I meant to ask about int(1/(sin(x)-1/2),x=Pi/4..Pi/2).
--------------------------------------------------------------------------------
> int(1/(sin(x)-1/2),x=Pi/4..Pi/2);
1/2 1/2 1/2 1/2 1/2
- 4/3 3 arctanh(1/3 3 ) - 4/3 3 arctanh(1/3 (2 - 3) 3 )
--------------------------------------------------------------------------------
> evalf(");
2.083882042
--------------------------------------------------------------------------------
