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