### Kristen Lew, Homework 12 ## ExpMath Spring 2012 Help:=proc(): print(`GuessDE(L,x,Di,Max), GuessDE1(L,x,Di,M), AlgToDe(F,P,x,Di,K), AlgToSeq(F,P,x,K,f0)`): end: ### Problem 3 ## GuessDE(L,x,Di,Max) that inputs a list of L, representing the first nops(L) terms of # a formal power series P(x)=L[1]+L[2]*x+..+L[nops(L)]x^(nops(L)-1)+.. # and outputs a differential operator with polynomial coefficients OPER(x,Di) # (where Di is differentiation with respect to x) such that the ORDER ("degree" in Di) and # DEGREE (degree in x) satisfy DEGREE + ORDER <= Max but with DEGREE+ORDER a small as # possible, such that when you apply OPER(x,Di) to P(x) it is (conjecturally) 0 (and surely # the first L coefficients of the Maclaurin expansion of OPER(x,Di) should all vanish. # Otherwise it is to return FAIL. # For example, GuessDE([seq(1/i!,i=0..20)],x,Di,1); # should return Di-1, since the e^x is a solution of the differential equation (Di-1)e^x=0. ### I couldn't figure this out... Di:=proc(F,n,x): if n=0 then F: else diff(F,x$n): fi: end: GuessDE1:=proc(L,x,Di,M) local P, OPER, ORDER, DEGREE, m, i, j, apply, var, var1,k,d,ak: P:=add(L[i]*x^(i-1), i=1..nops(L)): if M=0 then return P: fi: for ORDER from 0 to M-1 do DEGREE:=M-ORDER: OPER:=add(add(a[i,j]*x^i*Di^j, i=0..ORDER), j=0..DEGREE): var:={seq(seq(a[i,j], i=0..ORDER),j=0..DEGREE)}: apply:=convert(expand(OPER*P),list): for k from 1 to nops(apply) do if degree(apply[k],Di)>0 then d:=degree(apply[k],Di): ak:=apply[k]/Di^d: apply[k]:=Di(ak,degree(apply[k],Di),x): fi: od: apply:=expand(convert(apply,`+`)): var1:=solve(apply,var): apply:=expand(subs(var1, apply)): apply:=subs({seq(v=1, v in var)},apply): if apply=0 then OPER:=subs(var1,OPER): OPER:=subs({seq(v=1, v in var)},OPER): if OPER<>0 then RETURN(subs(var1,OPER)): fi: fi: od: FAIL: end: GuessDE:=proc(L,x,Di,Max) local k, g, P: P:=add(L[i]*x^(i-1), i=1..nops(L)): if Max=0 then return P: fi: for k from 1 to Max do g:=GuessDE1(L,x,Di,k): if g<>FAIL then RETURN(g): fi: od: FAIL: end: ### Problem 4 ## AlgToDe(F,P,x,Di,K) inputs a polynomial F(P,x) in the variables P and x, # representing an algebraic equation F(P(x),x)=0 satisfied by a formal power series P(x) # and a positive integer K and conjectures a differential operator by annihilating P(x) # by using the first K+1 coeffiencts of the Maclaurin expansion of P(x). # For example, AlgToDe(P*(1-x)-1,P,x,30); # should return (1-x)*Di-1 AlgToDe:=proc(F,P,x,Di,K) local L: L:=AlgToSeq(F,P,x,K,1): GuessDE(L,x,Di,K-1): end: AlgToSeq:=proc(F,P,x,K,f0) local H,F1,i,f,L,m: H:=coeff(F,P,1): if subs(x=0,H)=0 then RETURN(FAIL): fi: #a0(x)+H(x)*P(x)+a2(x)*P(x)^2+...=0 #P(x)=-(a0(x)+a2(x)*P(x)^2+....)/H(x) F1:=(P*H-F)/H: #start at f(x)=1 and iterate f(x)->F1(f(x)) f:=f0: L:=[f0]: for i from 1 to K do f:=subs(P=f,F1): f:=taylor(f,x=0,i+1): L:=[op(L),coeff(f,x,i)]: f:=add(coeff(f,x,m)*x^m,m=0..i): od: L: end: