csum:= proc(ep, np::name) # Test convergence of series # usage: # csum(ep, np) # ep an expression involving name np # determine convergence of sum(ep, np=a..infinity) # (for sufficiently large a that ep is defined for integers np >= a) # return true (converges), false (diverges), FAIL (undecided) # with infolevel[csum] >= 1, tell whether convergence is absolute # with infolevel[csum] >= 2, tell results of each test option `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; local e,n, ipars, k, limzero, eg1, eg2, ea, qpos, b; global `csum/a`, `property/object`; # assume n a large enough integer ipars:= map(x -> x^2, indets(ep, integer)); `csum/a`:= max(10000, op(ipars)); `property/object`[n]:= AndProp(RealRange(`csum/a`,infinity),integer); e:= eval(ep, np=n); userinfo(1, csum, `Checking sum of`,e,`over`,n); # sum of 0 converges absolutely; 0 would cause trouble in some tests if Testzero(e) then userinfo(2,csum, `Terms are identically 0`); userinfo(1,csum, `Converges absolutely`); return true fi; if hastype(e, negative^linear(n)) then # split up positive and negative `property/object`[k]:= AndProp(RealRange(`csum/a`,infinity),integer); eg1:= simplify(eval(e,n=2*k)); eg2:= simplify(eval(e,n=2*k+1)); if Testzero(eg1) or Testzero(eg2) then userinfo(3,csum,`Removing terms of 0`); return csum(eval(eg1+eg2),k); fi; userinfo(2,csum,`Alternate terms`,eg1,`and`,eg2); return `csum/alternate`(eg1,eg2,k); fi; # check limit of terms is 0 limzero:= `csum/limitzero`(e,n); if limzero = false then userinfo(1,csum,`Series diverges`); return false; fi; # check if terms are constant sign if is(e <= 0) then # if negative, change sign ea:= -e; qpos:= true; userinfo(3,csum,`Changing sign to make terms >= 0`); elif is(e >= 0) then ea:= e; qpos:= true; userinfo(3,csum,`Terms >= 0`); else ea:= `csum/abs`(e); qpos:= evalb(ea=e); if qpos then userinfo(3,csum,`Terms >= 0`); fi; fi; # Test absolute convergence b:= `csum/abscon`(ea,n); if b = true then userinfo(1,csum,`Series`,e,`converges absolutely`); return true; elif b = false then if qpos then userinfo(1,csum,`Series diverges`); return false; else userinfo(1,csum,`Series does not converge absolutely`); fi elif qpos then userinfo(1,csum,`Could not decide convergence`) else userinfo(1,csum,`Could not determine absolute convergence`); fi; # Test conditional convergence `csum/condit`(e,n); end; `csum/alternate`:= proc(e1, e2, n) # test convergence of series where terms are alternately e1 and e2 # usage: # `csum/alternate`(e1, e2, n) # e1 and e2: expressions in n # n: a name, assumed to be a sufficiently large positive integer # series is e1(n=a) + e2(n=a) + e1(n=a+1) + e2(n=a+1) + ... # return true (series converges), false (diverges), FAIL (undecided) option `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; local es, ea, limzero, i, qpos, qch, b; es:= [e1, e2]; # check limits of terms as n -> infinity # if either is nonzero, diverges limzero:= [FAIL,FAIL]; for i from 1 to 2 do limzero[i]:= `csum/limitzero`(es[i],n); if limzero[i] = false then userinfo(1,csum,`Series diverges`); return false; fi; od; # check if terms are constant sign ea:= [0,0]; qpos:= [FAIL,FAIL]; qch:= [false,false]; for i from 1 to 2 do if is(es[i] <= 0) then # if negative, change sign ea[i]:= -es[i]; qpos[i]:= true; qch[i]:= true; userinfo(3,csum,`Changing sign in`,es[i],`to make terms >= 0`); elif is(es[i] >= 0) then ea[i]:= es[i]; qpos[i]:= true; userinfo(3,csum,`Terms`,es[i],`>= 0`); else ea[i]:= `csum/abs`(es[i]); qpos[i]:= evalb(ea[i]=es[i]); if qpos[i] then userinfo(3,csum,`Terms`,es[i],`>= 0`); fi; fi; od; # Test absolute convergence # series converges absolutely iff both converge absolutely # if all terms have same sign, series diverges if either does not conv.abs. b:= map(`csum/abscon`,ea,n); if b = [true,true] then userinfo(1,csum,`Series converges absolutely`); return true; fi; if member(false,b) then if qpos = [true,true] and qch[1]=qch[2] then userinfo(1,csum,`Series diverges`); return false; else userinfo(1,csum,`Series does not converge absolutely`); fi else userinfo(1,csum,`Could not determine absolute convergence`); fi; # Test conditional convergence # if both converge, series converges # if one converges and other diverges, series diverges for i from 1 to 2 do if (b[i] = FAIL) and not qpos[i] then b[i]:= `csum/condit`(es[i],n); fi od; if b = [true,true] then userinfo(1,csum,`Series`,e,`converges`); return true; elif b = [true,false] or b = [false,true] then userinfo(1,csum,`Series diverges`); return false; fi; # test convergence of e1+e2 series if at least one limit known to be 0 if member(true,limzero) then userinfo(2,csum,`Regrouping`); csum(e1+e2,n); else userinfo(1,csum,`Could not determine convergence`); FAIL fi end; `csum/limitzero`:= proc(e,n) # test whether limit(e,n=infinity) = 0 # returns true, false or FAIL # note that a nonconstant result is assumed nonzero # even if it's 0 for some values of parameters # e: an expression in n # n: name, assumed positive integer option `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; local L, ar, br; userinfo(2, csum, `Checking limit of terms`); L := traperror(limit(e, n = infinity)); if L <> lasterror then if Testzero(L) then userinfo(2, csum, `OK, limit is 0`); return true fi; if (type(L, algebraic) and not has(L,'limit')) or (typematch(L, (ar::constant) .. (br::constant)) and is(0 < ar or br < 0)) then userinfo(1, csum, `Diverges, limit of terms is `, L); return false; fi; fi; L := traperror(limit(`csum/abs`(e), n = infinity)); if L <> lasterror then if Testzero(L) then userinfo(2, csum, `OK, limit is 0`); return true fi; if (type(L, algebraic) and not has(L,'limit')) or (typematch(L, (ar::constant) .. (br::constant)) and is(0 < ar or br < 0)) then userinfo(1, csum, `Diverges, limit of absolute values of terms is `, L); return false; fi; fi; userinfo(1,csum,`limit of absolute values of terms is`, L); FAIL end; `csum/abs`:= proc(x) local r, v, k, u; option `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; if is(0 <= x) then return x elif is(x <= 0) then return -x end if; r := abs(x); if has(r, abs) then if typematch(r, abs(v::anything)) then if type(v, `*`) then return map(`csum/abs`, v) elif typematch(v, (u::algebraic)^(k::algebraic)) and is(k, real) then return `csum/abs`(u)^k end if elif type(r, `*`) then return map(`csum/abs`, r) end if end if; r end proc; `csum/abscon`:= proc(e, n) # check absolute convergence of sum # usage: # `csum/abs`(e, n) # e an expression involving name n # where e >= 0 for sufficiently large integer n option `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; local b, test; for test in [`csum/ratiotest`, `csum/loglogtest`, `csum/indef`, `csum/split`] do b:= test(e,n); if b <> FAIL then return b fi od; FAIL end; `csum/condit`:= proc(e,n) # check for convergence, terms may not be nonnegative # this is not called if terms are known to be nonnegative # e: an expression depending on name n option `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; local test, b; for test in [`csum/osctest`,`csum/indef`, `csum/split`] do b := test(e, n); if b = true then userinfo(1, csum, `Series converges`); return true; elif b = false then userinfo(1, csum, `Series diverges`); return false; fi; od; userinfo(1, csum, `Could not determine convergence`); FAIL end proc; `csum/ratiotest`:= proc(e,n) # ratio test # e an expression in variable n, which should be nonnegative for large n local L, ep, b; option `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; userinfo(2, csum, `Trying ratio test`); L := traperror( limit(subs(n = n + 1, e)/e, n = infinity)); if L = lasterror or has(L, limit) then if typematch(e,abs(ep::algebraic)) then L:= traperror(limit(subs(n=n+1,ep)/ep, n=infinity)); else ep:= e fi; if L = lasterror or has(L, limit) then b:= simplify(subs(n=n+1,ep)/ep); L := traperror(limit(b, n = infinity)) fi fi; if type(L, range) then b := traperror(is(-1 < op(1, L)) and is(op(2, L) < 1)); if b = true then userinfo(2, csum, `Ratio test succeeds with ratio limit `, L); return true end if; b := traperror(is(1 < op(1, L)) or is(op(2, L) < -1)); if b = true then userinfo(2, csum, `Ratio test fails with ratio limit `, L); return false end if else # not a range if L = lasterror or L = 1 or L = -1 then userinfo(2, csum, `Ratio test inconclusive with ratio limit`, L); return FAIL end if; b := traperror(is(L < 1) and is(-1 < L)); if b = true then userinfo(2, csum, `Ratio test succeeds with ratio limit `, L); return true end if; b := traperror(is(1 < L) or is(L < -1)); if b = true then userinfo(2, csum, `Ratio test shows divergence with ratio limit`, L); return false end if end if; userinfo(2, csum, `Ratio test inconclusive with ratio limit`, L); FAIL end proc; `csum/loglogtest`:= proc(e,n) # loglog test # e an expression in variable n, which should be nonnegative for large n local L, b, lev, j, x, a; option `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; userinfo(2, csum, `Trying loglog test`); for lev to 10 do x:= ln(e*n*mul((ln@@j)(n), j = 1 .. lev - 1))/ (ln@@(lev + 1))(n); L := traperror(limit(x, n = infinity)); if L = lasterror or has(L,limit) then L:= traperror(limit(asympt(x,n), n=infinity)); fi; if L = lasterror or has(L,limit) then if is(x < -1) then L:= -2; userinfo(4,csum, ln(a[n]*n*mul((ln@@j)(n), j = 1 .. lev - 1))/ (ln@@(lev + 1))(n),`< -1 as `,n -> infinity); elif is(x > -1) then L:= 0; userinfo(4,csum, ln(a[n]*n*mul((ln@@j)(n), j = 1 .. lev - 1))/ (ln@@(lev + 1))(n),`> -1 as `,n -> infinity); fi else userinfo(4, csum, `limit of `, ln(a[n]*n*mul((ln@@j)(n), j = 1 .. lev - 1))/ (ln@@(lev + 1))(n), is, L); fi; if L <> -1 then break end if end do; if type(L, range) then b := traperror(is(op(2, L) < -1)); if b = true then userinfo(2, csum, `Loglog test succeeds`); return true end if; b := traperror(is(-1 < op(1, L))); if b = true then userinfo(2, csum, `Loglog test shows divergence`); return false end if elif L = -1 then userinfo(2, csum, `Loglog test fails`); return FAIL elif L <> lasterror then b := traperror(is(L < -1)); if b = true then userinfo(2, csum, `Loglog test succeeds`); return true end if; b := traperror(is(-1 < L)); if b = true then userinfo(2, csum, `Loglog test shows divergence`); return false end if end if; userinfo(2, csum, `Loglog test inconclusive`); FAIL end proc; `csum/indef`:= proc(e,n) # indefinite summation test: if you can find an indefinite # sum, see if it has a limit local F, L; option `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; userinfo(2, csum, `Trying indefinite summation`); F := sum(e, n); if has(F, 'sum') then userinfo(2, csum, `Indefinite summation failed`); return FAIL end if; L := limit(F, n = infinity); if has(L, 'limit') or type(L, range) then userinfo(2, csum, `Could not determine limit: failed`); return FAIL end if; if hastype(L, infinity) or type(L, undefined) then userinfo(2, csum, `Limit of indefinite sum is`, L); return false end if; userinfo(2, csum, `Indefinite summation succeeded`); true end proc; `csum/split`:= proc(e,n) # split up the summand into pieces and see if the pieces # converge local i, b, bd, bf, bs, qpos; option `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; if not type(e, `+`) then return FAIL end if; userinfo(2, csum, `Splitting up terms of`, e); b := convert(e, list); bs:= convert(map(signum,b),set); qpos:= evalb((bs = {1}) or (bs = {-1})); # all positive or all negative: if any diverges the sum does b := map(csum, b, n); bd := 0; bf := 0; for i to nops(b) do if b[i] = false then if qpos then userinfo(2, csum, `Term has divergent sum`); return false fi; bd := bd + 1 elif b[i] = FAIL then bf := bf + 1 end if end do; if bd = 1 and bf = 0 then userinfo(2, csum, `Sum of one term diverges`); return false end if; if bd = 0 and bf = 0 then userinfo(2, csum, `Sums of all terms converge`); true else userinfo(2, csum, `Splitting up failed`); FAIL end if end proc; `csum/osctest`:= proc(e, n, a) local ep, u, v, f, S, L, q; option `Maple Advisor Database 1.01 for Maple 6`, `Copyright (c) 2000 by Robert B. Israel. All rights reserved`; # oscillating series test # sum(u[n]*v[n],n) converges if # sum(u[n],n) is bounded and v[n] -> 0 monotonically ep := normal(e); if type(ep, `*`^algebraic) then ep := normal(expand(e)) end if; if not type(ep, `*`) then return FAIL end if; userinfo(2, csum, `Trying oscillating series test on `, ep); u := 1; v := 1; for f in ep do q:= indets(f, (-1)^algebraic); if q <> {} then u:= u*f; next fi; q:= indets(f, specfunc(anything, exp)); if (q <> {}) and traperror(is(I*op(op(1,q)),real)) = true then u:= u*f; next fi; q:= indets(f, specfunc(anything,{cos,sin})); if (q <> {}) and traperror(is(op(op(1,q)),real)) = true then u:= u*f; next fi; v := f*v; end do; userinfo(4, csum, Splitting, e, into, u, `and`, v); S := sum(u, n); if has(S, 'sum') then userinfo(4, csum, `Could not explicitly sum `, u); userinfo(2, csum, `Oscillating series test fails`); return FAIL end if; L := traperror(limit(S, n = infinity)); if type(L, infinity) then userinfo(4, csum, `Indefinite sum`, S, of, u, `has infinite limit`); userinfo(2, csum, `Oscillating series test fails`); return FAIL end if; if has(L, 'limit') or L = lasterror then userinfo(4, csum, `Could not find limit of indefinite sum`, S) ; userinfo(2, csum, `Oscillating series test fails`); return FAIL end if; userinfo(4, csum, `Indefinite sum`, S, `is bounded`); L := traperror(limit(v, n = infinity)); if L <> 0 then userinfo(4, csum, `Could not show convergence to 0 of`, v); userinfo(2, csum, `Oscillating series test fails`); return FAIL end if; L := traperror(asympt(v, n)); if L = lasterror then userinfo(4, csum, `Could not check monotonicity of`, v); userinfo(2, csum, `Oscillating series test fails`); return FAIL end if; L := traperror(limit(L/`csum/abs`(L), n = infinity)); if L = 1 or L = -1 then userinfo(2, csum, `Oscillating series test succeeds`); return true end if; userinfo(4, csum, `Could not show monotonicity of`, v); FAIL end proc;