#!/usr/local/bin/maple # -*- maplev -*- # Nathaniel Shar # HW 15 - 16 Mar 2014 # Experimental Mathematics # Stuff from c15.txt: #C15.txt: How Runge and Kutta could have (and should have #discovered (and proved!) their revolutionary method #March 23, 2014 Help15:=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 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 #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:=simplify(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 Help:=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)