#Nathan Fox #Homework 14 #I give permission for this work to be posted online #Read procedures from class read(`C14.txt`): Help:=proc(): print(` DelinZheng(lambda) , RP(x,y,M,d,N:=1) , srProveB(B,d,M,K,N:=1) `): end: ##PROBLEM 1## #convinced ##PROBLEM 2## #ProveB([ [[],[1/2],[0,1/2],[0,0,1]],[1/6,1/3,1/3,1/6],[0,1/2,1/2,1] ],4); returns true #ProveB([ [[],[1/2],[0,1/2],[0,0,1]],[1/4,1/4,1/4,1/4],[0,1/2,1/2,1] ],4); returns false, actually order 2 #ProveB([ [[],[1/2],[0,1/2],[0,0,1]],[1/10,2/10,3/10,4/10],[0,1/2,1/2,1] ],4); returns false, actually order 1 #ProveB([ [[],[1/2],[0,1/2],[0,0,1]],[1/10,4/10,4/10,1/10],[0,1/2,1/2,1] ],4); returns false, actually order 2 ##PROBLEM 3## #DelinZheng(lambda): inputs a number lambda, and outputs the #Runge-Kutta method described at the bottom of p.1 and the top of #p.2 of the wikipedia article, as a triplet [a,b,c], the so-called #Butcher matrix. DelinZhang:=proc(lambda) return [[[], [1/2], [1/2-1/lambda, 1/lambda], [0, 1-lambda/2, lambda/2]], [1/6, (4-lambda)/6, lambda/6, 1/6], [0, 1/2, 1/2, 1]]: end: ##PROBLEM 4## #ProveB(DelinZhang(i), 4); returned true for all integers #i from 1 to 5 #ProveB(DelinZhang(2.5), 4); returned false, but this was due to #roundoff error, because ProveB(DelinZhang(5/2), 4); returned #true. This method has order 4. #The same story holds for 3.5 and 7/2 #In fact, ProveB(DelinZhang(lambda), 4); returns true, so this is #a proof for all lambda ##PROBLEM 5## #RP(x,y,M,d,N:=1): inputs variables x,y, a positive integer M, and #a positive integer d and outputs a polynomial in x,y, of (total) #degree with numeric, integer coefficients between N and M #(note: N is optional and defaults to 1, but it's important to #be able to test with negative coefficients sometimes) RP:=proc(x, y, M, d, N:=1) local i,j,r: r:=rand(N..M): return add(add(r()*x^i*y^j, j=0..d-i), i=0..d): end: ##PROBLEM 6## #srProveB(B,d,M,K,N:=1): rigorously proves that B is NOT of order #d if it returns false, and semi-rigorously proves that B is of #order d if it returns true, by checking K different randomly #chosen polynomials generated by RP(x,y,M,d,N). #(Note: N is optional and defaults to 1, but just in case you want #to test negative coefficients, you can do so) srProveB:=proc(B,d,M,K,N:=1) local f,h,a,Y1,Y2,Y,i: for i from 1 to K do f:=RP(x,y,M,d,N): Y1:=GRK1(op(B),f,x,y,0,0,h): Y2:=subs(x=h,Z(f,x,y,0,d)): Y:=expand(Y1-Y2): if ldegree(Y,h)