Physics 350:
Computational Methods
for the Physical Sciences

Spring Semester 2009

Solutions for Physics 350 Lab 2:
Taylor's Series and Complex Arithmetic
 

As per usual, we will start by clearing Maple's memory: 

> restart;

Part II: Taylor's Series in Maple

First of all we are asked to calculate the Taylor's series of some functions as a way of verifying Euler's formula (at least around
t = 0).  We'll do this by doing Taylor's series expansions to 15th order of exp(I*t), cos(t), and sin(t) and showing exp(I*t) = cos(t) + I*sin(t). 

> f := t -> exp(I*t);  # one way to create a function of the variable t
proc (t) options operator, arrow; exp(`*`(I, `*`(t))) end proc (1.1)
 
> expt := taylor(f(t),t=0,15);   # 15th order Taylor's series expansion about t=0
series(`+`(1, `*`(I, `*`(t)), `-`(`*`(`/`(1, 2), `*`(`^`(t, 2)))), `-`(`*`(`*`(`/`(1, 6), `*`(I)), `*`(`^`(t, 3)))), `*`(`/`(1, 24), `*`(`^`(t, 4))), `*`(`*`(`/`(1, 120), `*`(I)), `*`(`^`(t, 5))), `-`...
series(`+`(1, `*`(I, `*`(t)), `-`(`*`(`/`(1, 2), `*`(`^`(t, 2)))), `-`(`*`(`*`(`/`(1, 6), `*`(I)), `*`(`^`(t, 3)))), `*`(`/`(1, 24), `*`(`^`(t, 4))), `*`(`*`(`/`(1, 120), `*`(I)), `*`(`^`(t, 5))), `-`...
series(`+`(1, `*`(I, `*`(t)), `-`(`*`(`/`(1, 2), `*`(`^`(t, 2)))), `-`(`*`(`*`(`/`(1, 6), `*`(I)), `*`(`^`(t, 3)))), `*`(`/`(1, 24), `*`(`^`(t, 4))), `*`(`*`(`/`(1, 120), `*`(I)), `*`(`^`(t, 5))), `-`...
(1.2)
 
> expt2 := convert(expt,polynom);   # Converting the results into a polynomial
`+`(1, `*`(I, `*`(t)), `-`(`*`(`/`(1, 2), `*`(`^`(t, 2)))), `-`(`*`(`+`(`*`(`/`(1, 6), `*`(I))), `*`(`^`(t, 3)))), `*`(`/`(1, 24), `*`(`^`(t, 4))), `*`(`*`(`/`(1, 120), `*`(I)), `*`(`^`(t, 5))), `-`(`...
`+`(1, `*`(I, `*`(t)), `-`(`*`(`/`(1, 2), `*`(`^`(t, 2)))), `-`(`*`(`+`(`*`(`/`(1, 6), `*`(I))), `*`(`^`(t, 3)))), `*`(`/`(1, 24), `*`(`^`(t, 4))), `*`(`*`(`/`(1, 120), `*`(I)), `*`(`^`(t, 5))), `-`(`...
`+`(1, `*`(I, `*`(t)), `-`(`*`(`/`(1, 2), `*`(`^`(t, 2)))), `-`(`*`(`+`(`*`(`/`(1, 6), `*`(I))), `*`(`^`(t, 3)))), `*`(`/`(1, 24), `*`(`^`(t, 4))), `*`(`*`(`/`(1, 120), `*`(I)), `*`(`^`(t, 5))), `-`(`...
(1.3)
 
> g := t -> cos(t);
proc (t) options operator, arrow; cos(t) end proc (1.4)
 
> cost := convert(taylor(g(t),t=0,15),polynom);
`+`(1, `-`(`*`(`/`(1, 2), `*`(`^`(t, 2)))), `*`(`/`(1, 24), `*`(`^`(t, 4))), `-`(`*`(`/`(1, 720), `*`(`^`(t, 6)))), `*`(`/`(1, 40320), `*`(`^`(t, 8))), `-`(`*`(`/`(1, 3628800), `*`(`^`(t, 10)))), `*`(... (1.5)
 
> temp := sin(t); h := unapply(temp,t);  # creating a function of t another way
sin(t)
 

 

proc (t) options operator, arrow; sin(t) end proc (1.6)
 
> sint := convert(taylor(h(t),t=0,15),polynom);    
`+`(t, `-`(`*`(`/`(1, 6), `*`(`^`(t, 3)))), `*`(`/`(1, 120), `*`(`^`(t, 5))), `-`(`*`(`/`(1, 5040), `*`(`^`(t, 7)))), `*`(`/`(1, 362880), `*`(`^`(t, 9))), `-`(`*`(`/`(1, 39916800), `*`(`^`(t, 11)))), ... (1.7)
 
> euler_test := expt2 - (cost + I*sint);
`+`(`*`(`*`(`/`(1, 362880), `*`(I)), `*`(`^`(t, 9))), `-`(`*`(`+`(`*`(`/`(1, 39916800), `*`(I))), `*`(`^`(t, 11)))), `*`(I, `*`(t)), `*`(`*`(`/`(1, 6227020800), `*`(I)), `*`(`^`(t, 13))), `-`(`*`(`+`(...
`+`(`*`(`*`(`/`(1, 362880), `*`(I)), `*`(`^`(t, 9))), `-`(`*`(`+`(`*`(`/`(1, 39916800), `*`(I))), `*`(`^`(t, 11)))), `*`(I, `*`(t)), `*`(`*`(`/`(1, 6227020800), `*`(I)), `*`(`^`(t, 13))), `-`(`*`(`+`(...
(1.8)

Notice that euler_test doesn't look like zero... but all the terms cancel out.  There are a couple of ways to get this reduced to zero.  You can "simplify" the expersion: 

> simplify(euler_test);
0 (1.9)

Or you could evaluate the complex number using evalc: 

> evalc(euler_test);
0 (1.10)

Or you can convert euler_test into a function of t and collect the terms: 

> euler_test_func := unapply (euler_test,t);
proc (t) options operator, arrow; `+`(`*`(`*`(`/`(1, 362880), `*`(I)), `*`(`^`(t, 9))), `-`(`*`(`+`(`*`(`/`(1, 39916800), `*`(I))), `*`(`^`(t, 11)))), `*`(I, `*`(t)), `*`(`*`(`/`(1, 6227020800), `*`(I...
proc (t) options operator, arrow; `+`(`*`(`*`(`/`(1, 362880), `*`(I)), `*`(`^`(t, 9))), `-`(`*`(`+`(`*`(`/`(1, 39916800), `*`(I))), `*`(`^`(t, 11)))), `*`(I, `*`(t)), `*`(`*`(`/`(1, 6227020800), `*`(I...
(1.11)
 
> collect(euler_test_func(t),t);
0 (1.12)
 
And so we see that at least to 15th order, the Taylor's series expansions around zero of exp(I*t), cos(t), and sin(i) are compliant with Euler's formula. 

Next we are asked to calculate the Taylor's series of cos(x) about x=0 and compare it (in a plot) to cos(x) in the range x=0 to 5.
 
> j := x -> cos(x);
proc (x) options operator, arrow; cos(x) end proc (2.1)
 
> jtay := convert(taylor(j(x),x=0,10),polynom);
`+`(1, `-`(`*`(`/`(1, 2), `*`(`^`(x, 2)))), `*`(`/`(1, 24), `*`(`^`(x, 4))), `-`(`*`(`/`(1, 720), `*`(`^`(x, 6)))), `*`(`/`(1, 40320), `*`(`^`(x, 8)))) (2.2)
 
> k := unapply(jtay,x);
proc (x) options operator, arrow; `+`(1, `-`(`*`(`/`(1, 2), `*`(`^`(x, 2)))), `*`(`/`(1, 24), `*`(`^`(x, 4))), `-`(`*`(`/`(1, 720), `*`(`^`(x, 6)))), `*`(`/`(1, 40320), `*`(`^`(x, 8)))) end proc (2.3)
 
> plot([j(x),k(x)],x=0..5,color=[red,blue]);# Plot the two functions together to see how they compare
Plot_2d  

Now we are asked to check at what x values the disagreement between these two functions is 0.1.  The suggestion is to use the fsolve function instead of solve since solve will just revert to fsolve when an algebraic solution is not found. 

> fsolve(k(x)-j(x)=0.1,x);
3.632624500 (3.1)
 

Now we are done with Taylor's series expansions, let's move onto complex arithmetic. 


Part II: Complex Number Arithmetic in Maple

As part 2 of this lab exercise, we are asked to perform some complex arithmetic in Maple.  The first thing I am asked to do is solve for the three roots of (1 + I)^(1/3).  It is suggested that we use the solve command rather than direct evaluation, since direct evaluation will only give on of the three roots. 

> m := (1+I)^(1/3);
`*`(`^`(`+`(1, I), `/`(1, 3))) (4.1)

> simplify(evalc(m));  # First we do this the one way we are told not to, a direct evaluation.  This only finds one root.
`*`(`^`(2, `/`(1, 6)), `*`(`+`(cos(`+`(`*`(`/`(1, 12), `*`(Pi)))), `*`(I, `*`(sin(`+`(`*`(`/`(1, 12), `*`(Pi))))))))) (4.2)

> evalf(%);
`+`(1.084215081, `*`(.2905145554, `*`(I))) (4.3)
 
> solve(x^3=(1+I),x);  # A better way to solve this...
`*`(`^`(`+`(1, I), `/`(1, 3))), `+`(`-`(`*`(`/`(1, 2), `*`(`^`(`+`(1, I), `/`(1, 3))))), `*`(`*`(`/`(1, 2), `*`(I)), `*`(`^`(3, `/`(1, 2)), `*`(`^`(`+`(1, I), `/`(1, 3)))))), `+`(`-`(`*`(`/`(1, 2), `*... (4.4)
 
> evalf(%);
`+`(1.084215081, `*`(.2905145555, `*`(I))), `+`(`-`(.7937005258), `*`(.7937005257, `*`(I))), `+`(`-`(.2905145552), `-`(`*`(1.084215081, `*`(I)))) (4.5)
 

 Now I am asked to solve for (4+I)^(4-I), which the writeup notes is a messy result: 

> n := (4+I)^(4-I);
`^`(`+`(4, I), `+`(4, `-`(I))) (5.1)
 
> p := simplify(evalc(n));
`+`(`-`(`*`(289, `*`(exp(arctan(`/`(1, 4))), `*`(`+`(`-`(cos(`+`(`*`(`/`(1, 2), `*`(ln(17))), `-`(`*`(4, `*`(arctan(`/`(1, 4)))))))), `*`(I, `*`(sin(`+`(`*`(`/`(1, 2), `*`(ln(17))), `-`(`*`(4, `*`(arc... (5.2)
 

This is just one of the roots, trying to use solve instead: 

> n_alt := solve(x^(1/(4-I)) = (4+I),x);
exp(`*`(`+`(4, `-`(I)), `*`(ln(`+`(4, I))))) (5.3)
 
> p_alt := simplify(evalc(n_alt));
`+`(`-`(`*`(289, `*`(exp(arctan(`/`(1, 4))), `*`(`+`(`-`(cos(`+`(`*`(`/`(1, 2), `*`(ln(17))), `-`(`*`(4, `*`(arctan(`/`(1, 4)))))))), `*`(I, `*`(sin(`+`(`*`(`/`(1, 2), `*`(ln(17))), `-`(`*`(4, `*`(arc... (5.4)
 Notice that the directly evaluated "p" is the same as the solved "p_alt", suggesting there is indeed only one solution.  Just to complete this, I solved for the real and imaginary portions of (4+I)^(4-I) and evaluated them as "floats" so see what x+I*y actually looked like. 

> re_p:=Re(p); evalf(re_p);
 
`+`(`*`(289, `*`(exp(arctan(`/`(1, 4))), `*`(cos(`+`(`*`(`/`(1, 2), `*`(ln(17))), `-`(`*`(4, `*`(arctan(`/`(1, 4))))))))))) 
334.5750532 (5.5)
 
> im_p := Im(p); evalf(im_p);
`+`(`-`(`*`(289, `*`(exp(arctan(`/`(1, 4))), `*`(sin(`+`(`*`(`/`(1, 2), `*`(ln(17))), `-`(`*`(4, `*`(arctan(`/`(1, 4))))))))))))  
-156.1614520 (5.6)
 
Now as a final task, I am asked to compute the arcsine of 16 and then confirm that when I add 2*Pi and take the sine that the result behaves as expected, despite the complex number nature of the result.  Let's do it... 
> q := arcsin(16);
arcsin(16) (6.1)
 
> simplify(evalc(q));
`+`(`*`(`/`(1, 2), `*`(Pi)), `-`(`*`(`+`(I), `*`(ln(`+`(16, `*`(`^`(255, `/`(1, 2))))))))) (6.2)
 
> r := q + 2*Pi;
`+`(arcsin(16), `*`(2, `*`(Pi))) (6.3)
 
> s := sin(r);
16 (6.4)