#hw12.txt; March 3, 2014; Frank Wagner Help:=proc(): print(` Zg1(f,x,y,cur) `): print(` Zg(f,x,y,y0,N `): print(` ERC(f,x) `): print(` PolAppx(f,x,y,y0,x1,N,eps) `): end: ################### #####Problem 1##### ################### Zg1:=proc(f,x,y,cur,n) local cur1,c: if n=ERC(Zg(f,x,y,y0,N)[N+1]/2,x)[1] then RETURN(FAIL): fi: while abs(subs(x=x1,Zg(f,x,y,y0,N+i)[N+i+1]-Zg(f,x,y,y0,N+i)[N+i]))>=eps/100 do i:=i+1: od: P:=evalf(Zg(f,x,y,y0,N+i)[N+i+1]): v:=subs(x=x1,P): RETURN(P,v): end: #PolAppx returns both the polynomial and its value at x=x1 since I'm not sure #if we wanted both or just the value. # # #PolAppx(1+x^2+y^3,x,y,1,0.3,50,10^(-10)); returns an incredibly long polynomial of #degree 116, whose value at x1 is: # 2.603375587 # #PolAppx(cos(x)*exp(y)+1+x+y,x,y,1,0.3,50,10^(-10)); returns # FAIL #since ERC(Zg(cos(x)*exp(y)+1+x+y,x,y,1,50)[51]/2,x) returns # [.2622802598, .2623901394], #and abs(x1)=0.3 is larger than both of those. ################### #####Problem 3##### ################### #Since MacLaurin series are polynomials, it's clear that any solution can be written as # Sum(a(n)*x^n,n=0..infinity). #It is obvious that the coefficient of x^0 is 1 by the initial condition. #So, a(n)=1=(1/0!) for n=0. Hence, this satisfies the form we desire. #Now, suppose a(n)=(c(n)/n!), where c(n) is an integer for all n less than or equal to some positive integer k. #Then, as we saw in problem 1, the coefficient a(k+1) of k+1 in the MacLaurin polynomial #will fulfill (a(k+1))*(k+1)=coeff(taylor(subs(y=cur1,f),x=0,k+1),x,k). #By the way Taylor/MacLaurin polynomials are built, we can see that the right side is # [F^(k)](0)/(k!), where F=subs(y=cur1,f) and F^(k) denotes the kth derivative of F. #Multiplying the k! to each side, we get (a(k+1))*((k+1)!)=[F^(k)](0) #Since a(k+1) is already paired with x^(k+1) and we're dealing with a polynomial, the power of x that is #multiplying by a(k+1) in F(x) can only increase. Hence, when we take the kth derivative of x, #a(k+1) will still have a power of x multiplying by it, which will vanish when we evaluate F^(k) at 0. #Now, because the coefficient for every power of x is an integer in f and in cur (and clearly the powers #themselves are integers), we can see that [F^(k)](0) is some integer. Call it c(k+1). #So, a(k+1)=(c(k+1))/((k+1)!), fitting our desired form. #By induction, we're done. ################### #####Problem 4##### ################### Zm1:=proc(f,x,y,cur) local n,cur1,c,d,p,q,k,i,j,fn,g,var,eq: if nops(cur)<>nops(f) then RETURN(FAIL): fi: if nops(cur)<>nops(y) then RETURN(FAIL): fi: n:=[]: for i from 1 to nops(cur) do if cur[i]=0 then d:=0: else d:=degree(cur[i],x): fi: n:=[op(n),d]: od: c:=[]: for i from 1 to nops(cur) do q:=k[i]: c:=[op(c),q]: od: cur1:=[]: for i from 1 to nops(cur) do p:=cur[i]+c[i]*x^(n[i]+1): cur1:=[op(cur1),p]: od: fn:=[]: for i from 1 to nops(cur) do g:=subs({seq(y[j]=cur1[j],j=1..nops(cur))},f[i]): fn:=[op(fn),g]: od: var:={seq(c[i],i=1..nops(c))}: eq:={seq( coeff(diff(cur1[i],x)-subs({seq(y[j]=cur1[j],j=1..nops(cur1))},f[i]),x,n[i])=0,i=1..nops(cur1))}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: for i from 1 to nops(cur) do c[i]:=eval(c[i],var[i]): od: for i from 1 to nops(cur) do cur1[i]:=cur[i]+c[i]*x^(n[i]+1): od: cur1: end: Zm:=proc(f,x,y,y0,N) local i,L,l,cur,p,n,cur1: L:=[]: cur:=[]: for i from 1 to nops(f) do p:=y0[i]: cur:=[op(cur),p]: od: for i from 1 to nops(f) do l:=y0[i]: cur1:=cur: for n from 1 to N do cur1:=Zm1(f,x,y,cur1): l:=[op(l),cur1[i]]: od: L:=[op(L),l]: od: L: end: #Zm([1+x*y1+y2, 3+x*y1^2+y2],x,[y1,y2],[1,2],10); returned: # [[1, 1+3*x, 3*x^2+3*x+1, 2*x^3+3*x^2+3*x+1, 2*x^3+3*x^2+3*x+1+(3/2)*x^4, # 2*x^3+3*x^2+3*x+1+(3/2)*x^4+(13/10)*x^5, 2*x^3+3*x^2+3*x+1+(3/2)*x^4+(13/10)*x^5+(17/15)*x^6, # 2*x^3+3*x^2+3*x+1+(3/2)*x^4+(13/10)*x^5+(17/15)*x^6+(53/60)*x^7, # 2*x^3+3*x^2+3*x+1+(3/2)*x^4+(13/10)*x^5+(17/15)*x^6+(53/60)*x^7+(437/672)*x^8, # 2*x^3+3*x^2+3*x+1+(3/2)*x^4+(13/10)*x^5+(17/15)*x^6+(53/60)*x^7+(437/672)*x^8+(2873/6048)*x^9, # 2*x^3+3*x^2+3*x+1+(3/2)*x^4+(13/10)*x^5+(17/15)*x^6+(53/60)*x^7+(437/672)*x^8+(2873/6048)*x^9+(53107/151200)*x^10], # [2, 2+5*x, 3*x^2+5*x+2, 3*x^3+3*x^2+5*x+2, 3*x^3+3*x^2+5*x+2+(9/2)*x^4, 3*x^3+3*x^2+5*x+2+(9/2)*x^4+(53/10)*x^5, # 3*x^3+3*x^2+5*x+2+(9/2)*x^4+(53/10)*x^5+(293/60)*x^6, 3*x^3+3*x^2+5*x+2+(9/2)*x^4+(53/10)*x^5+(293/60)*x^6+(1709/420)*x^7, # 3*x^3+3*x^2+5*x+2+(9/2)*x^4+(53/10)*x^5+(293/60)*x^6+(1709/420)*x^7+(3799/1120)*x^8, # 3*x^3+3*x^2+5*x+2+(9/2)*x^4+(53/10)*x^5+(293/60)*x^6+(1709/420)*x^7+(3799/1120)*x^8+(86549/30240)*x^9, # 3*x^3+3*x^2+5*x+2+(9/2)*x^4+(53/10)*x^5+(293/60)*x^6+(1709/420)*x^7+(3799/1120)*x^8+(86549/30240)*x^9+(717071/302400)*x^10]] # #Zm([x+y1,y1*y2+x,x+y1*y2*y3],x,[y1,y2,y3],[1,2,1],10); returned: # [[1, 1+x, x^2+x+1, x^2+x+1+(1/3)*x^3, x^2+x+1+(1/3)*x^3+(1/12)*x^4, x^2+x+1+(1/3)*x^3+(1/12)*x^4+(1/60)*x^5, # x^2+x+1+(1/3)*x^3+(1/12)*x^4+(1/60)*x^5+(1/360)*x^6, x^2+x+1+(1/3)*x^3+(1/12)*x^4+(1/60)*x^5+(1/360)*x^6+(1/2520)*x^7, # x^2+x+1+(1/3)*x^3+(1/12)*x^4+(1/60)*x^5+(1/360)*x^6+(1/2520)*x^7+(1/20160)*x^8, # x^2+x+1+(1/3)*x^3+(1/12)*x^4+(1/60)*x^5+(1/360)*x^6+(1/2520)*x^7+(1/20160)*x^8+(1/181440)*x^9, # x^2+x+1+(1/3)*x^3+(1/12)*x^4+(1/60)*x^5+(1/360)*x^6+(1/2520)*x^7+(1/20160)*x^8+(1/181440)*x^9+(1/1814400)*x^10], # [2, 2+2*x, 2+2*x+(5/2)*x^2, 2+2*x+(5/2)*x^2+(13/6)*x^3, 2+2*x+(5/2)*x^2+(13/6)*x^3+(11/6)*x^4, # 2+2*x+(5/2)*x^2+(13/6)*x^3+(11/6)*x^4+(22/15)*x^5, 2+2*x+(5/2)*x^2+(13/6)*x^3+(11/6)*x^4+(22/15)*x^5+(13/12)*x^6, # 2+2*x+(5/2)*x^2+(13/6)*x^3+(11/6)*x^4+(22/15)*x^5+(13/12)*x^6+(1927/2520)*x^7, # 2+2*x+(5/2)*x^2+(13/6)*x^3+(11/6)*x^4+(22/15)*x^5+(13/12)*x^6+(1927/2520)*x^7+(10469/20160)*x^8, # 2+2*x+(5/2)*x^2+(13/6)*x^3+(11/6)*x^4+(22/15)*x^5+(13/12)*x^6+(1927/2520)*x^7+(10469/20160)*x^8+(61547/181440)*x^9, # 2+2*x+(5/2)*x^2+(13/6)*x^3+(11/6)*x^4+(22/15)*x^5+(13/12)*x^6+(1927/2520)*x^7+(10469/20160)*x^8+(61547/181440)*x^9+(97261/453600)*x^10], # [1, 1+2*x, 1+2*x+(9/2)*x^2, 1+2*x+(9/2)*x^2+(47/6)*x^3, 1+2*x+(9/2)*x^2+(47/6)*x^3+(27/2)*x^4, # 1+2*x+(9/2)*x^2+(47/6)*x^3+(27/2)*x^4+(263/12)*x^5, 1+2*x+(9/2)*x^2+(47/6)*x^3+(27/2)*x^4+(263/12)*x^5+(2435/72)*x^6, # 1+2*x+(9/2)*x^2+(47/6)*x^3+(27/2)*x^4+(263/12)*x^5+(2435/72)*x^6+(126667/2520)*x^7, # 1+2*x+(9/2)*x^2+(47/6)*x^3+(27/2)*x^4+(263/12)*x^5+(2435/72)*x^6+(126667/2520)*x^7+(729313/10080)*x^8, # 1+2*x+(9/2)*x^2+(47/6)*x^3+(27/2)*x^4+(263/12)*x^5+(2435/72)*x^6+(126667/2520)*x^7+(729313/10080)*x^8+(291737/2880)*x^9, # 1+2*x+(9/2)*x^2+(47/6)*x^3+(27/2)*x^4+(263/12)*x^5+(2435/72)*x^6+(126667/2520)*x^7+(729313/10080)*x^8+(291737/2880)*x^9+(62781949/453600)*x^10]]