# Solutions to Assignment 5 # # 1. Since x exp(-x^2) has an elementary antiderivative, we write # the integrand as u dv with u=1/x and dv = x exp(-x^2) dx. > with(student): > J:= exp(R^2) * Int(exp(-x^2), x=R .. inf > inity); > J1:= intparts(J, 1/x); # Note that for large R, the integral term here will be smaller than # J because of the 1/x^2. This time, take u = 1/x^3. > J1:= intparts(J1,1/x^3); # Now a pattern is starting to emerge: each time we get another x^2 # in the denominator. > for jj from 3 to 5 do J1:= intparts(J1, > 1/x^(2*jj-1)) od: > J1; # (b) After n iterations, the integral term is c(n) exp(-x^2)/x^(2 n), # where c(n) is a constant. Integrating that by parts: > intparts(c(n)*Int(exp(-x^2)/x^(2*n),x=R. > .infinity), 1/x^(2*n+1)); # So c(n+1) = -(2n+1)/2 c(n). Since c(0) = 1, we can write > c:= n -> (-1)^n*2^(-n) * Product(2*k+1, > k = 1 .. n-1); # Maple actually can find a "closed form" for this product. > value(c(n)); > c1:= simplify("); # Alternatively, note that c(n) contains the product of the odd # integers from 1 to 2n-1. The product of all the integers from 1 to # 2n is (2n)!, while the product of the even ones is # 2^n n!. So we can also write c(n) as > c2:= (-1)^n*2^(-2*n)*(2*n)!/n!; > simplify(c2/c1); > c:= unapply(c2,n); > seq(c(nn), nn = 1 .. 5); # Whichever way we write it, the n'th term of the series is > nthterm:=unapply( c(n-1)/2 * R^(-2*n+1), > n); > theseries:= Sum(nthterm(n), n=1..infinit > y); # Check this for the first 5 terms. > partialsum:= n -> sum(nthterm(k), k=1..n > ); > partialsum(5); # To check for convergence, use the Ratio Test. > nthterm(n+1)/nthterm(n); > simplify("); # For any R the absolute value of this goes to infinity as n -> infinity, # so the series diverges. # # (c) The remainder after 5 terms is > R5:= expand(J1 - partialsum(5)); # Another integration by parts yields > intparts(R5,x^(-11)); > expand("); # Since the integrand is positive, we conclude that # 0 > R5 > -945/64*R^(-11). For R=2 we have > evalf(945/64 * 2^(-11)); # We calculated nthterm(n+1)/nthterm(n) = - (2 n - 1)/(2 R^2). So # |nthterm(n+1)| >= |nthterm(n)| when 2 n - 1 >= 2 R^2. For R = 2 # this means n >= 5. > seq(evalf(abs(subs(R=2,nthterm(nn)))), n > n=1..6); # The remainder is smaller in absolute value than the next term of the # series. This estimate for the error is smallest for the fourth partial # sum. # Let's check the actual errors. > V:= value(J); > seq(evalf(subs(R=2,V - partialsum(nn))), > nn = 1 .. 5); # # 2. We'll try it with the partial sum to n=10, and up to 10 terms of # the Euler-Maclaurin series. > f:= n -> 1/(n^2+1): > PartialSum:= N -> Sum(f(n), n = 1 .. N): > readlib(eulermac): Digits:= 20: > p:= evalf(PartialSum(10)); > ems:= [seq(eulermac(f(N),N,2*dd-1), dd= > 1..10)]: > ems; # Note that this has a limit of Pi/2 as N -> infinity. So if the sum of # the series is S, F(N) = PartialSum(N-1) - S + Pi/2. Our # approximation # for S will be PartialSum(10) - F(11) + Pi/2). > ems:= evalf(eval(subs(N=11,O=0,ems))); # S to S all agree to within 10^(-9). The best value might be # S, since it and S agree to the accuracy of our calculation. # So our answer would be > p - ems+evalf(Pi/2); # Compare to Maple's value for the sum. > sum(f(n), n=1..infinity); > evalf("); # 3. (a) > read `fib.def`: Warning: new definition for norm Warning: new definition for trace mpow := proc(n) options remember; if type(n,posint) then if type(n,even) then evalm(mpow(1/2*n)^2) else evalm(mpow(1/2*n-1/2)^2 &* M ) fi elif type(n,negint) then evalm(1/mpow(-n)) elif type(n,`+`) then evalm(`&*`(op(map(mpow,n)))) elif type(n,`*`) then evalm(mpow(n/op(1,n))^op(1,n)) else matrix([ [F[n-1],F[n]],[F[n],F[n]+F[n-1]] ]) fi end zeck := proc(x) local n; if x = 0 then 0 else n := invfib(x); F[n]+zeck(x-fib(n)) fi end > fib(n+4)+fib(n-4); # So the equation is F[n+4] + F[n-4] = 7 F[n]. # (b) The identity was > fib(m+n); # We conclude that F[m+n] is divisible by F[n] if F[m] is. Since # F = 0 is divisible by F[n], we conclude by mathematical # induction that F[k n] is divisible by F[n] for all positive integers k. # # (c) Since F = 3, all the Fibonacci numbers F[4 k] are divisible # by 3. What's the first one divisible by 9? > seq(fib(4*kk)/9, kk = 1 .. 10); # So F is divisible by 9, and then F[12 k] is divisible by 9 for all # positive integers k. Which are divisible by 3^3 = 27? > seq(fib(12*kk)/27, kk = 1 .. 10); # F is the first, and then of course F[36 k] for all k. # 4, 12, 36: is there a pattern here? It looks like F[4*3^m] will be # divisible by 3^(m+1). If F[n] is divisible by 3^m, is F[3*n] # divisible by 3^(m+1)? > fib(3*n); > subs(F[n] = 3^m * r, "); > "/3^(m+1); > simplify("); # Yes: if m >= 1 then this is an integer. We conclude that # F[4*3^m*k] is divisible by 3^(m+1) for all nonnegative integers m # and k. # # (d) If F[n]/F[m] = 7880997 and n and m are not too small, we # can approximate F[n] and F[m] by fapp(n) and fapp(m) to get > fapp(n)/fapp(m); # So 7880997 should be approximately ((1+sqrt(5))/2)^(n-m). # Thus n - m should be approximately > evalf(ln(7880997)/ln((1+sqrt(5))/2)); # It looks very much like n-m = 33. If so, . > shouldbe:= fib(m+33)/fib(m)= 7880997; > solve(", F[m]); # So F[m] should be 89/55*F[m-1]. In particular, F[m-1] should be # divisible by 55. Is 55 itself a Fibonacci number? > invfib(55); > fib(10); > fib(11); # So we want m=11, making F[m] = 89 and F[m-1]=55, and # n=11+33=44. > fib(44)/fib(11); # # 4. We try a Chebyshev series. > with(numapprox): Digits:= 15: > f: x -> ln(1+x^3): > ff:= x -> evalf(ln(1+x^3)): > chebyshev(ff,x=0 .. 3, 10^(-5)); > ch:= unapply(eval(subs(T=orthopoly[T],") > ),x); > plot(ln(1+x^3)-ch(x), x = 0 .. 3); # According to the graph, the maximum error is approximately # 3.4*10^(-6) , which is acceptable. # Note: you could expand ch(x) into a sum of coefficients times # powers of x, but this would # be bad for numerical accuracy when x is near 3, so it's better not # to.