#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: 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: