General RK Scheme Order

Suppose we have the ODE in (ℝ^n)
(ODE)	y'[t]==f[t,y[t]]
and we use the s-stage RK scheme described by the Butcher tableau
	({{(Subscript[c, i]), (Subscript[A, i,j])}, { , (Subscript[b, j])}})
which means
(RK)	({{(Subscript[k, i]), =, f[(Subscript[t, in])+(Subscript[c, i])h,(Subscript[y, in])+h(∑_(j=1)^s)(Subscript[A, i,j])(Subscript[k, j])]}, {(Subscript[y, out]), =, (Subscript[y, in])+h (∑_(j=1)^s)(Subscript[b, j])(Subscript[k, j])}})
and can be evaluated directly in the event that the s×s array is strictly lower triangular.

Note:  If A has non-zero diagonal entries or upper triangular entries the k-equations are a nonlinear system for s vectors of size n.

All we need to do to get the equations that the coefficients A, b, and c need to satisfy for RK to have a specified order is form the power series for the LTE
	LTE=y[(Subscript[t, out])]-(Subscript[y, out])=y[(Subscript[t, in])+h]-(Subscript[y, out])
Then collect the powers of h and make sure that all the terms less than our desired order are identically zero for all functions.

There is a lot more book keeping here than you might think and there is a historical mistake that persisted for a long time because matrix products do not commute.  The derivation was originaly carried out for scalar equations where all the products commute and the results applied to vector equations where the products did not.

The power series computation for y[(Subscript[t, in])+h] is relatively straightforward.  Write down the series then differentiate the ODE to replace the higher derivatives of (y^((n)))[(Subscript[t, in])]. the following is to convince you that we do not want to do this computation for the non-autonomous system

Clear[Sub,n,y,f,t]
Sub[1]=y'[t]->f[t,y[t]];
Sub[n_]:=Sub[n]=Simplify[D[Sub[n-1],t]/.Sub[1]];
Sub[4]/.{
f[t,y[t]]->f,
Derivative[p1_,p2_][f][t,y[t]]->Derivative[p1,p2][f]
}

So we are going to focus on the autonomous system

Clear[Sub,n,y,f,t]
Sub[1]=y'[t]->f[y[t]];
Sub[n_]:=Sub[n]=D[Sub[n-1],t]/.Sub[1];
Sub[5]/.{
f[y[t]]->f,
Derivative[p1_][f][y[t]]->Derivative[p1][f]
}

This would be fine BUT Mathematica is assuming that products commute (we could turn this off but it is probably better to think about what is going one) and is tidying up the expressions for us.  If we used these we would be about to recreate a historical mistake!  Anyway for now take a look at the pieces of these expressions and realize they are all different and could be turned on/off for different functions f. A generic ODE solver has to deal with all functions and so it has to eliminate all of these independent pieces.

What Butcher and others did was look at these expressions in the messy non-autonomous case and and work out how you get them.  They built this strange looking construct on top of a tree structure that lists of all the independent pieces that you get for a general RK scheme and gives an easy way to write down the equations. We are going to look at this briefly for the non-autonomous system.

Clear[Sub,n,y,f,t]
Tidy={
f[t,y[t]]->f,
Derivative[p1_,p2_][f][t,y[t]]->Derivative[p1,p2][f]
};
Sub[1]={y'[t]->f[t,y[t]]};
Sub[n_]:=Sub[n]=D[Sub[n-1],t]/.Sub[1];

For the non-autonomous system the y'' term has two independent pieces and the order should be  
	Derivative[1,0][f]+Derivative[0,1][f]f 
where the gradient Derivative[0,1][f] is a matrix and f is a vector. See below.

Sub[2]/.Tidy

The y''' term is the result of differentiating the y'' term 
	Derivative[2,0][f]+Derivative[1,1][f]f+Derivative[1,1][f]f f+Derivative[0,2][f]f f + Derivative[0,1][f]Derivative[0,1][f] f   
where the second gradient Derivative[0,3][f] is a three indexed object.

The y'''' term is the result of differentiating the y''' term 
	 
where the second gradient Derivative[0,3][f] is a three indexed object.

Sub[3]/.Tidy