#!/usr/local/bin/maple # -*- maplev -*- # Nathaniel Shar # HW 14 - 12 Mar 2014 # Experimental Mathematics # 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 #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,x,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: ############# # Problem 2 # ############# # I always got order 1 (except for 1/6, 1/3, 1/3, 1/6, of course!) ############# # Problem 3 # ############# 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: # Yes, it's true for all values of lambda (even irrational ones like # Pi). You won't find this out if you set lambda to a floating-point number, # though! ############# # Problem 4 # ############# RP := proc(x, y, M, d) local r: if not type(M, integer) then: error("Type error"): fi: r := rand(1..M): add(add(r()*x^i*y^j,j=0..d-i),i=0..d): end: ############# # Problem 5 # ############# srProveB := proc(B, d, M, K) local i, P, x, y, Y, h: for i from 1 to K do: P := RP(x, y, M, d): Y := GRK1s(op(B), P, x, y, 0, 0, h, d) - subs(x=h, Z(P, x, y, 0, d)): if ldegree(Y, h) < d+1 then: print(sprintf("Conjecture disproved by polynomial %s", convert(P, string))): return false: fi: od: return true: end: