#hw15.txt, 3/15/14, Anthony Zaleski Help:=proc(): HelpC15(): print(`FindRK(s,d,a,b,constraints)`): end: ##########Problem 1: Finding general RK of order d########## #NOTE: I modified GButcher so that the b's sum to 1. #The method here uses 5 random polynomials to guess #a possible formula for an s-stage degree d RK method. #We must still check the formula with ProveB once it's #generated. #NOTE: var names a and b are also inputted since there may #be free variables in the solution, so it wouldn't do to make a,b local. #Additionally, one can optionally include constraints in #a,b to limit the number of solutions outputted. (See comments below.) FindRK:=proc(s,d,a,b,constraints:={}) local f,B,h,i,j,var: #f:=RP(x,y,3,d): B:=GButcher(a,b,s): var:={seq(seq(a[i,j],j=1..i-1),i=2..s)} union {seq(b[i],i=1..s-1)}: subs(solve({seq(seq(coeff(DiffRK(B,d,RP(x,y,3,d),x,y,h),h,i),i=0..d),j=1..5)} union constraints,var),B): end: #I was able to use this to find the general form of an #order 2 method with 2 stages (which I verified with ProveB): #F:=FindRK(2,2,a,b); #F := [[[], [a[2, 1]]], [(1/2)*(2*a[2, 1]-1)/a[2, 1], 1-(1/2)*(2*a[2, 1]-1)/a[2, 1]], [0, a[2, 1]]] #As you can see, the formula is alread complicated. The #general formula for a (3,3) method is even worse; it involves #roots which Maple can't find. I ran out of time trying a general #(4,4) method. #However, if I add the constraints {a[3,1]=0,a[3,2]=1/2} to the equations #in the solve part and look for (4,4) methods, I get RK4: #F:=FindRK(4,4,a,b,{a[3,1]=0,a[3,2]=1/2}); #F := [[[], [1/2], [0, 1/2], [0, 0, 1]], [1/6, 1/3, 1/3, 1/6], [0, 1/2, 1/2, 1]] #With two different constraints, we can recover a Delin-Zheng method: #F:=FindRK(4,4,a,b,{a[3,1]=1/6,a[3,2]=1/3}); #F := [[[], [1/2], [1/6, 1/3], [0, -1/2, 3/2]], [1/6, 1/6, 1/2, 1/6], [0, 1/2, 1/2, 1]] #However, it takes too long to recover ALL D-Z methods by passing a #symbolic constraint in terms of lambda. #So I guess my procedure is only useful if we have sufficient constraints #on the tableau; otherwise there are way too many possibilities. #For s=4, it seems like specifying a[3,2] and a[3,3] is usually enough. #The problem is, just using arbitrary constraints doesn't always produce #an interesting result. Since it helps to know the actual method before #hand to give the right constraints, this kind of defeats the purpose... #Well, at least we found a formula for second-order R-K methods! #In fact, after the simplification below, we see that we've #derived Wikipedia's formula for second-order R-K! #see http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#Second-order_methods_with_two_stages #subs(a[2,1]=alpha,simplify(FindRK(2,2,a,b))); #[[[], [alpha]], [(1/2)*(2*alpha-1)/alpha, (1/2)/alpha], [0, alpha]] ############################################################### #################C15.txt WITH GBUTCHER MODIFIED SO b's SUM TO 1 ############################################################### #C15.txt: How Runge and Kutta could have (and should have #discovered (and proved!) their revolutionary method #March 23, 2014 HelpC15:=proc(): print(` DelinZheng(lambda) `): print(`srProveB(B,d,M,K,N:=1)`): print(`GButcher(a,b,s) , DiffRK(B,d,f,x,y,h) `): end: #stuff from C14.txt #C14.txt, March 10, 2014, Runge-Kutta methods HelpC14:=proc(): print(`GRK1(a,b,c,f,x,y,x0,y0,h)`): print(`GRK(a,b,c,f,x,y,x0,y0,h,x1)`): print(`Z(f,x,y,y0,N), P(x,y,a,d) `): print(`ProveB([ [[],[1/2],[0,1/2],[0,0,1]],[1/6,1/3,1/3,1/6],[0,1/2,1/2,1] ],4);`): end: RK4:= [[],[1/2],[0,1/2],[0,0,1]],[1/6,1/3,1/3,1/6],[0,1/2,1/2,1]: #GRK1(a,b,c,f,x,y,x0,y0,h): One step in the #General Runge-Kutta with Butcher matrix a,b,c #(see wiki) GRK1:=proc(a,b,c,f,x,y,x0,y0,h) local k,s,i: s:=nops(b): for i from 1 to s do k[i]:=h*subs({x=x0+c[i]*h,y=y0+add(a[i][j]*k[j],j=1..i-1)},f): od: y0+add(b[i]*k[i],i=1..s): end: #GRK(a,b,c,f,x,y,x0,y0,h,x1): The general Runge-Kutta #with Butcher matrix a,b,c appx. to #the ode #y'(x)=f(x,y(x)), y(x0)=y0 all the way #to x1 with step-size h #the output is pair of lists #Xvalues, Yavlues GRK:=proc(a,b,c,f,x,y,x0,y0,h,x1) local xc,yc,Lx,Ly,s,i: if not (type(a,list) and type(b,list) and type(c,list)) then RETURN(FAIL): fi: s:=nops(b): if nops(a)<>s or nops(c)<>s then RETURN(FAIL): fi: if add(b[i],i=1..s)<>1 then RETURN(FAIL): fi: if c[1]<>0 then RETURN(FAIL): fi: Lx:=[]: Ly:=[]: xc:=x0: yc:=y0: while xc<=x1 do Lx:=[op(Lx),xc]: Ly:=[op(Ly),yc]: yc:=GRK1(a,b,c,f,x,y,xc,yc,h): xc:=xc+h: od: Lx,Ly: end: ##from C12.txt #Z1(f,x,y,cur): if cur is the current Maclaurin #polynomial approximating y'(x)=f(x,y(x)) finds #the poly appx. of one degree higher #in the Picard sequence that is supposed to #converge to the solution of y'(x)=f(x,y(x)), y(x0)=y0 #Right now f is a polynomial of x and y #Version of March 7, 2014 correcting semantical bugs spotted by #Andrew Lohr Z1:=proc(f,x,y,cur) local n,cur1,c: if cur=0 then n:=0: else n:=degree(cur,x): fi: cur1:=cur+c*x^(n+1): c:=solve(coeff(diff(cur1,x)-subs(y=cur1 ,f),x,n)=0,c): cur+c*x^(n+1): end: #Z(f,x,y,y0,N): the first N+1 MacLaurin polynomials #n x of sol. of y'(x)=f(x,y(x)), y(0)=y0 #Version of March 7, 2014 correcting semnatical bugs spotted by #Andrew Lohr Z:=proc(f,x,y,y0,N) local cur,n: if diff(f,y)=0 then print(`the function f only depends on x, so the solution is the indefinite integral`): RETURN(y0+int(f,x=0..x)): fi: cur:=y0: for n from 1 to N do cur:=Z1(f,x,y,cur): od: cur: end: ###End of stuff from C12.txt #P(x,y,a,d): a generic polynomial in x,y of degree d #in variables x,y and degree d P:=proc(x,y,a,d) local i,j: add(add(a[i,j]*x^i*y^j,j=0..d-i),i=0..d): end: #GRK1s(a,b,c,f,x,y,x0,y0,h,R): One step in the #General Runge-Kutta with Butcher matrix a,b,c #to order R, done symbolically (h is symbolic!) #and at every step we truncate above h^R #(see wiki) GRK1s:=proc(a,b,c,f,x,y,x0,y0,h,R) local k,s,i,j: s:=nops(b): for i from 1 to s do k[i]:=expand(h*subs({x=x0+c[i]*h,y=y0+add(a[i][j]*k[j],j=1..i-1)},f)): k[i]:=expand(add(coeff(k[i],h,j)*h^j,j=0..R)): od: expand(y0+add(b[i]*k[i],i=1..s)): end: #ProveB(B,d): inputs a Butcher triplet [a,b,c] #and a positive integer d, rigorously proves that the #explicit Runge-Kutta method, based on B is of order d, #try: #ProveB([ [[],[1/2],[0,1/2],[0,0,1]],[1/6,1/3,1/3,1/6],[0,1/2,1/2,1] ],4); ProveB:=proc(B,d) local f,h,a,Y1,Y2,Y: f:=P(x,y,a,d): Y1:=GRK1s(op(B),f,x,y,0,0,h,d): Y2:=subs(x=h,Z(f,x,y,0,d)): Y:=expand(Y1-Y2): if ldegree(Y,h)>=d+1 then true: else false: fi: end: ###end of stuff from C14.txt #C14.txt, March 10, 2014, Runge-Kutta methods Help14:=proc(): print(`GRK1(a,b,c,f,x,y,x0,y0,h)`): print(`GRK(a,b,c,f,x,y,x0,y0,h,x1)`): print(`Z(f,x,y,y0,N), P(x,y,a,d) `): print(`ProveB([ [[],[1/2],[0,1/2],[0,0,1]],[1/6,1/3,1/3,1/6],[0,1/2,1/2,1] ],4);`): end: RK4:= [[],[1/2],[0,1/2],[0,0,1]],[1/6,1/3,1/3,1/6],[0,1/2,1/2,1]: #GRK1(a,b,c,f,x,y,x0,y0,h): One step in the #General Runge-Kutta with Butcher matrix a,b,c #(see wiki) GRK1:=proc(a,b,c,f,x,y,x0,y0,h) local k,s,i: s:=nops(b): for i from 1 to s do k[i]:=h*subs({x=x0+c[i]*h,y=y0+add(a[i][j]*k[j],j=1..i-1)},f): od: y0+add(b[i]*k[i],i=1..s): end: #GRK(a,b,c,f,x,y,x0,y0,h,x1): The general Runge-Kutta #with Butcher matrix a,b,c appx. to #the ode #y'(x)=f(x,y(x)), y(x0)=y0 all the way #to x1 with step-size h #the output is pair of lists #Xvalues, Yavlues GRK:=proc(a,b,c,f,x,y,x0,y0,h,x1) local xc,yc,Lx,Ly,s,i: if not (type(a,list) and type(b,list) and type(c,list)) then RETURN(FAIL): fi: s:=nops(b): if nops(a)<>s or nops(c)<>s then RETURN(FAIL): fi: if add(b[i],i=1..s)<>1 then RETURN(FAIL): fi: if c[1]<>0 then RETURN(FAIL): fi: Lx:=[]: Ly:=[]: xc:=x0: yc:=y0: while xc<=x1 do Lx:=[op(Lx),xc]: Ly:=[op(Ly),yc]: yc:=GRK1(a,b,c,f,x,y,xc,yc,h): xc:=xc+h: od: Lx,Ly: end: ##from C12.txt #Z1(f,x,y,cur): if cur is the current Maclaurin #polynomial approximating y'(x)=f(x,y(x)) finds #the poly appx. of one degree higher #in the Picard sequence that is supposed to #converge to the solution of y'(x)=f(x,y(x)), y(x0)=y0 #Right now f is a polynomial of x and y #Version of March 7, 2014 correcting semantical bugs spotted by #Andrew Lohr Z1:=proc(f,x,y,cur) local n,cur1,c: if cur=0 then n:=0: else n:=degree(cur,x): fi: cur1:=cur+c*x^(n+1): c:=solve(coeff(diff(cur1,x)-subs(y=cur1 ,f),x,n)=0,c): cur+c*x^(n+1): end: #Z(f,x,y,y0,N): the first N+1 MacLaurin polynomials #n x of sol. of y'(x)=f(x,y(x)), y(0)=y0 #Version of March 7, 2014 correcting semnatical bugs spotted by #Andrew Lohr Z:=proc(f,x,y,y0,N) local cur,n: if diff(f,y)=0 then print(`the function f only depends on x, so the solution is the indefinite integral`): RETURN(y0+int(f,x=0..x)): fi: cur:=y0: for n from 1 to N do cur:=Z1(f,x,y,cur): od: cur: end: ###End of stuff from C12.txt #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. DelinZheng:=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: #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: #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)