#!/usr/local/bin/maple # -*- maplev -*- # Nathaniel Shar # HW 9 - 23 Feb 2014 # Experimental Mathematics with(LinearAlgebra): ############# # Problem 1 # ############# applyUmbra := proc(P, x, U): return add(U(i)*coeff(P, x, i), i=0..degree(P, x)): end: ############# # Problem 2 # ############# Pn := proc(n, x, a, b) local P, A, eqns, vars, solns: option remember: P := add(A[i]*x^i, i=0..n): eqns := {seq(int(P*x^i, x=a..b)=0, i=0..n-1), A[n]=1}: vars := {seq(A[i], i=0..n)}: solns := solve(eqns, vars): return subs(solns, P): end: PnG := proc(n, x, U) local P, A, eqns, vars, solns: option remember: P := add(A[i]*x^i, i=0..n): eqns := {seq(applyUmbra(P*x^i, x, U)=0, i=0..n-1), A[n]=1}: vars := {seq(A[i], i=0..n)}: solns := solve(eqns, vars): return subs(solns, P): end: ############# # Problem 3 # ############# GuessRFd := proc(L, dtop, dbot, n) local R, N, D, i, eqns, vars, solns: if nops(L) <= dtop + dbot +7 then: error "Not enough terms to guess": fi: N := add(A[i]*n^i, i=0..dtop): D := add(B[i]*n^i, i=0..dbot): eqns := {seq(denom(L[i])*subs({n=i}, N) = numer(L[i])*subs({n=i}, D), i=1..nops(L)), A[dtop] = 1}: vars := {seq(A[i], i=0..dtop), seq(B[i], i=0..dbot)}: solns := solve(eqns, vars): if solns = NULL then: return FAIL: fi: R := subs(solns, N/D): if {seq(simplify(subs({n=i}, R)- L[i]), i=1..nops(L))} <> {0} then: return FAIL: fi: return factor(normal(R)): end: GuessRF := proc(L, n) local cand, dtop, dbot: dtop := 0: dbot := 0: while true do: cand := GuessRFd(L, dtop, dbot, n, r): if cand <> FAIL then: return cand: else: if dtop = 0 then: dtop := dbot + 1: dbot := 0: else: dtop := dtop - 1: dbot := dbot + 1: fi: fi: od: end: GuessCoeff := proc(n, i, U, terms := 20) local x, L: L := [seq(coeff(PnG(j, x, U), x, j-i), j=1..terms)]: return GuessRF(L, n): end: # See next part for general answers. ############# # Problem 4 # ############# GuessDoubleRF := proc(L, m, n): return GuessRF([seq(GuessRF(L[i], n), i=1..nops(L))], m): end: Ratios := proc(L) local i: return [seq(L[i+1]/L[i], i=1..nops(L)-1)]: end: GuessOP := proc(U, terms := 20) local n, r, x, L, i, t, base, rat, rat2, fac, success: base := GuessCoeff(n, 0, U): L := Ratios([seq(GuessCoeff(n, i, U, 2*terms+8), i=0..terms)]): rat := GuessRF(L, t): # Unfortunately, Maple is VERY STUPID when it comes to # multiplication, and does not know how to do # even extremely simple things like prod((n-t), t=1..r) when n-t # is a positive integer, which is pathetic, so in order to get a # closed form I would have to write a new simplification engine # for products of rational functions. # I am not about to do that, so we'll just return a function # instead of an expression. Everything will work nicely as long as # you don't try to plug in symbolic values for the second argument # (which is "r"). If you do that, you might get something in terms # of gamma functions with negative integer arguments. If that # happens, then you could use the (nonsensical) simplification # rule # GAMMA(b-a)/(GAMMA(-a)*GAMMA(b+1)) = (-1)^b*b!*binomial(a, b) # if a-b is a nonnegative integer. return (x, y)->subs({n=x}, base) * product(subs({n=x, t=j}, rat), j=1..y): end: