High Order Explicit RK Schemes

Richard Fears and Jason Gregersen

Purpose

To explore the construction of high order explicit RK schemes

Introduction

The modern organization of RK schemes presented by Butcher ([1] pg 94,95) utilizes a Butcher tableau in the following form. (({{c, A}, {0, b}})) This leads to the following algorithm for an RK scheme. (Subscript[y, n+1])=(Subscript[y, n])+h (∑_(i=1)^s)(Subscript[b, i]) (Subscript[k, i]) where: (Subscript[k, i])=f(yn+h(∑_(j=1)^s)(Subscript[a, ij])(Subscript[k, j]))

Thus to determine an RK scheme one just needs to construct the proper entries into the table to be used in the algorithm.  The number of rows of A is the number of “stages” of the scheme. The order of the scheme is related to the order of accuracy for the scheme.  In general we need at least as many stages as the order, but as the desired order increases the number of stages required becomes larger than the desired order. Coefficients that lead to schemes of a given order are found by satisfying a set of equations knows as order conditions.

Order Conditions

To establish the correct order conditions we can compare the Taylor series expansion of the exact solution and the approximate solution. It is through this process that theory supporting the appropriate order conditions was established. As a method to actually calculating the appropriate order conditions we can look at the “Butcher trees”. Each order is associated with a family of trees. The tree has as many vertices as the order. Each tree leads to an equation. Thus to form the order conditions for a third order RK scheme one has to look at the family of trees with one vertex and form the appropriate order condition, then move on to the family of trees with two vertices and the family of trees with three vertices. The process of creating an order condition from a specific tree is detailed at the following website; http://www.123mathvids.com/vids/rkordereqn/rkordereqn.html  As an example, the order conditions for a third order scheme with “s” stages is shown below:

<<NumericalDifferentialEquationAnalysis`
Eqs=RungeKuttaOrderConditions[3,s];

Eqs

Here we can see that only one condition was formed for the first and second order, as there exists only one tree in each of those groups. Using this process we can construct a set of order equations, but since this becomes a large set of nonlinear equations quickly (for a 12th order method there are 7813 order conditions [3]) it becomes difficult to solve this system and then construct the associated high order RK schemes. Additional conditions are generally added, but an optimal strategy for all cases is not defined.

Project Outline

1. We will construct and solve the order equations needed to form a third order scheme. 

2. We will discuss reasons for choosing choosing parameters to form specific schemes and compare our choice to those of others.

3.  We will perform the same analysis for a 5th order 6 stage scheme.

4.  We will discuss trying to solve larger systems by forming a minimization problem.

Code

Function

f[{y1_,y2_,y3_,y4_}]:=
Module[{α},
α=((((y1^2)+(y2^2)))^(3/2));
{y3,y4,-(y1/α),((-y2)/α)}
];

Initial Values

xfinal=40.*π//N;
steps=20000;
e=.75;
posy0={1-e,0,0,(Sqrt[((1+e)/(1-e))])};
vinit={-1-e,0,0,-(Sqrt[((1-e)/(1+e))])};
h=0.01;
Tol=1.*(10^(-10));

Tableau's from the previous homework

Euler={(({{0, 0}, {1, 0}})),(({{1, 0}, {1, -1}}))};
Heun={(({{0, 0}, {1, 0}})),(({{(1/2), (1/2)}, {1, -1}}))};

Bogacki={{{0,0,0,0},{1/2,0,0,0},{0,3/4,0,0},{2/9,1/3,4/9,0}},{{2/9,1/3,4/9,0},{7/24,1/4,1/3,1/8}}};

Fehlberg={{{0,0,0,0,0,0},{1/4,0,0,0,0,0},{3/32,9/32,0,0,0,0},{1932/2197,−7200/2197 ,7296/2197,0,0,0},{439/216,−8, 3680/513,−845/4104,0,0},{-8/27, 2,−3544/2565,1859/4104,−11/40,0}},{{25/216,0,1408/2565,2197/4104,−1/5,0},{16/135,0,6656/12825,28561/56430,−9/50,2/55}}};

DormandPrince={{{0,0,0,0,0,0,0},{1/5,0,0,0,0,0,0},{3/40,9/40,0,0,0,0,0},{44/45,-56/15,32/9,0,0,0,0},{19372/6561,-25360/2187,64448/6561,-212/729,0,0,0}, {9017/3168,-355/33,46732/5247,49/176,-5103/18656,0,0},{35/384,0,500/1113,125/192,-2187/6784,11/84,0}},{{5179/57600,0,7571/16695,393/640,-92097/339200,187/2100,1/40},{35/384,0,500/1113,125/192,-2187/6784,11/84,0}}};

CashKarp={{{0,0,0,0,0,0},{1/5,0,0,0,0,0},{3/40,9/40,0,0,0,0},{3/10,−9/10 ,6/5,0,0,0},{-11/54,5/2, -70/27,35/27,0,0},{1631/55296, 175/512,575/13824,44275/110592,253/4096,0}},{{37/378,0,250/621,125/594,0,512/1771},{2825/27648,0,18575/48384,13525/55296,277/14336,1/4}}};

General ODE solver

ODESolver is a function that takes in the tableau and function information and then takes a step of lentgh h from the point y. In my stepper i have created a matrix of k values "kvector". The dot product of this with the matrix from the tableau sets up my set of equations. and then i am making a table that sets the specific k values. The output is the higher order scheme and an error estimate.

ODESolver[Tableau_,y_,h_,f_]:=
Module[
{Schemesize,ynew,znew,kvector},
Schemesize=Length[Tableau[[1]]];
kvector=ConstantArray[0,Schemesize];
kvector[[1]]=f[y];Do[kvector[[i]]=f[y+h*Sum[kvector[[j]]*Tableau[[1]][[i,j]],{j,1,Schemesize}]],{i,2,Schemesize}];
ynew=y+h*Sum[Tableau[[2]][[1]][[i]]*kvector[[i]],{i,1,Schemesize}];
znew=y+h*Sum[Tableau[[2]][[2]][[i]]*kvector[[i]],{i,1,Schemesize}];
{ynew,Norm[ynew-znew]}
]

Step Size Adjuster

My step size adjuster below is taking in the error, a preset tolerance "Tol", and the order of the scheme. It then is computing the value "r" which is the ratio of the tolerance to the error. It then multiplies h by error and has some safety feature set so that i can't do more than double or half the step size in one step.

varstepsize[err_,Tol_,order_,hOld_]:=(
Module[{r,h},
r=((((Tol/err)))^(1/order));
r=Min[Max[.5,.9*r],2.0];
h=hOld*r;
h
])

Drivers

The driver below is doing the following:
1. Running the ODESolver and storing the results as ynew and err.
2. Says if the step is good take it and keep track of some stuff, and if it is bad don't take it and keep track of some stuff.
3. running the step through the step size program and then keeping track of the how far we have traveled with dist
4. Throught the loop i am sowing the data i want. The reap is then stored and data, and the final ouput is ...

varstepdriver[yinit_,Tableau_,Tol_,hinit_,f_]:=
Module[{hi=hinit,ye=yinit,dist=0,order=Length[Tableau[[1]]],stepstaken=0,totloops=0,k=0,ynew,err,data},
data=Reap[
While[dist<xfinal,
totloops++;
{ynew,err}=ODESolver[Tableau,ye,hi,f];
If[
err<Tol,
{Sow[ye=ynew,point],Sow[{dist,(((err))/(hi^order))},distvsdiff],Sow[{dist,err},distvserr],dist=dist+hi,stepstaken++},
{ye}];
hi=varstepsize[err,Tol,order,hi];
If[(*this step is just making sure the last step is the right size*)
(dist+hi)≥xfinal,
hi=(xfinal-dist)];
]];
(*{data,stepstaken,totloops}*)
{data[[2]][[1]][[All,1;;2]],data[[2]][[2]],data[[2]][[3]],stepstaken,totloops}
]

Below is the way the driver was structured for previous hw assignment. More information from the data was needed specifically the velocity information so I re-strucured the output of the driver.

Fixedstepdriver[Tableau_,yn_,f_,numsteps_,xfin_]:=
Module[{y,data,h},
h=(xfin/numsteps);
y=yn;
data=Table[y=ODESolver[Tableau,y,h,f][[1]],{numsteps}]
]

data=Fixedstepdriver[DormandPrince,posy0,f,1000,π];
ListPlot[data[[All,1;;2]]]

data=varstepdriver[posy0,DormandPrince,Tol,h,f];

data[[4]]

ListPlot[data[[1]]]

Misc Code

<<NumericalDifferentialEquationAnalysis`
OrderEqs[order_,stages_]:=Module[
{Eqs,rowsum},
Eqs=Flatten[RungeKuttaOrderConditions[order,stages,
ButcherSimplify->True,
RungeKuttaMethod->Explicit]];
rowsum=Table[Sum[ (Subscript[\[FormalA], i,j]),{j,1,i-1}]==(Subscript[\[FormalC], i]),{i,1,stages}];
Eqs=Join[Eqs,rowsum];
{Eqs,Length[Eqs],((((stages^2)-stages))/2)+2stages}
]
OrderEqsRow[order_,stages_]:=Module[
{Eqs,rowsum},
Eqs=Flatten[RungeKuttaOrderConditions[order,stages,
ButcherRowSum->True,
RungeKuttaMethod->Explicit]];
rowsum=Table[Sum[ (Subscript[\[FormalA], i,j]),{j,1,i-1}]==(Subscript[\[FormalC], i]),{i,1,stages}];
Eqs=Join[Eqs,rowsum];
{Eqs,Length[Eqs],((((stages^2)-stages))/2)+2stages}
]

Perr[order_,stages_,sol_]:=Module[
{err},
err=ButcherPrincipalError[order,stages,
ButcherSimplify->True,
RungeKuttaMethod->Explicit];
err=err/.sol
]

FormVar[order_,stages_,numvar_]:=Module[
{A,b,c},
A=LowerTriangularize[Table[(Subscript[\[FormalA], i,j]),{i,1,stages},{j,1,stages}],-1];
b=Table[(Subscript[\[FormalB], i]),{i,1,stages}];
c=Table[(Subscript[\[FormalC], i]),{i,1,stages}];
Table[ArrayRules[Flatten[{A,b,c}]][[i]][[2]],{i,1,numvar}]
]

FormTableau[order_,stages_,sol_]:=Module[
{A,b,c,b2},
b2=ConstantArray[0,stages];
A=LowerTriangularize[Table[(Subscript[\[FormalA], i,j]),{i,1,stages},{j,1,stages}],-1];
b=Table[(Subscript[\[FormalB], i]),{i,1,stages}];
{A/.sol,{b/.sol,b2}}
]

truesol[y1_]=y2/.Solve[{(((y1+e))^2)+((y2^2)/(1-(e^2)))==1},y2][[2]];//Quiet
AsymptoticRKErrorTest[y_,f_, Scheme_,Schemename_,{hMin_,hMax_},Color_] := Module[
{n,nMax=12,Data,h, yOut, Err,g,logh,appy1,appy2},
Data=Table[
h=hMin*((((hMax/hMin)))^(n/nMax));
{appy1,appy2}=ODESolver[Scheme,posy0,h,f][[1,1;;2]];
Err=Abs[truesol[appy1]-appy2];
{h,Err},
{n,0,nMax}
];
g[logh_]=Fit[Log[Data],{1, logh},logh];
ListLogLogPlot[Data,
PlotLabel->StringForm["``: Err = ``",Schemename,(E^(g[Log["h"]]))],
FrameLabel->{"h","err"},
Frame->True,PlotStyle->Color
]
]

(*To use the stability results here we are using the code from [2]*)
StabilityRegion[Scheme_,Color_]:=
Module[{g,λ,R,y0,z,y},
g[y_]=λ y;
R[z_]=Simplify[((ODESolver[Scheme,y0,(z/λ),g][[1]])/y0)];
ContourPlot[Norm[R[x+I y]]==1,{x,-6,4},{y,-4,4},PlotLabel->"Stability Region",ContourStyle->Color]
]
StabilityRegionEquation[Scheme_]:=
Module[{g,λ,R,y0,z,y},
g[y_]=λ y;
R[z_]=Simplify[((ODESolver[Scheme,y0,(z/λ),g][[1]])/y0)]
]

Quadsol[order_,stages_]:=
Module[{Vand,repl,γ,bvec,sys,bsol},
Vand=Array[((Subscript[\[FormalC], #2])^(#1-1))&,{order,stages}];
repl=Table[(Subscript[\[FormalC], i])->(i/stages),{i,2,stages}];
Vand=Vand/.repl;
Vand=Vand/.(Subscript[\[FormalC], 1])->0;
γ=Table[(1/(k+1)),{k,0,order-1}];
bvec=Table[(Subscript[\[FormalB], i]),{i,1,stages}];
sys=Table[(Vand.bvec)[[i]]==γ[[i]],{i,1,order}];
bsol=Solve[sys][[1]];
Join[bsol,repl]
]

Small Scheme Creation and Special Properties

Creating Small Schemes

In order to create an RK scheme we first need to establish and then solve the necessary order equations as derived earlier. We started with a straight forward third order three stage scheme. Adding in the standard Butcher Simplifications (column sum, row sum) and specifying that the scheme be explicit results in a nonlinear system of 8 equations and 9 variables. Given the small size of this system we are able to solve it directly using the “Solve” command in mathematica. The result is:

{order,stages}={3,3};
{Eqs,numeqs,numvar}=OrderEqs[order,stages];
sol=Solve[Eqs][[1]]

Choosing Parameters, error vs stability

Thus any choice of (Subscript[\[FormalC], 2]) should yield a third order RK scheme. We are now going to explore some options in choosing a value of (Subscript[\[FormalC], 2]). It is standard practice to calculate the Principal Error terms for this scheme. That is, to determine the accuracy of the scheme we look at the difference between one step of the scheme to the taylor series expansion at that point. Since our values satisfy the order equations for a third order scheme all of the terms to that power will be eliminated and the values in the next highest power (the first group not eliminated) are the principal error terms. When we take this error and replace the values with our solution we will get the principal err as a function of our parameter.

err=Perr[3,3,sol]

Error minimization

One method for choosing a value for (Subscript[\[FormalC], 2]) is to choose it so that it minimizes the error term. Typically this is done with the L2 Norm, but we would like choose a few others and compare the results. Using the L1,L2,L3, and L∞ we find the results for the value of the optimal (Subscript[\[FormalC], 2]) in each case. Since the difference between the L1, L2, and L3 norm are minimal. We will only compare the results using the L2 Norm, L∞ norm, and a case where we choose a value at random. We then put these optimal values back into our solution and form the appropriate Tableaus. Lastly we compare the results. By comparing our answer to the analytic solution to the Keplar two body problem for varying stepsizes. The results are show below.

smallerroptions={Minimize[Norm[err,1],(Subscript[\[FormalC], 2])][[2]],Minimize[Norm[err,2],(Subscript[\[FormalC], 2])][[2]],Minimize[Norm[err,3],(Subscript[\[FormalC], 2])][[2]],NMinimize[Norm[err,∞],(Subscript[\[FormalC], 2]),AccuracyGoal->20,WorkingPrecision->20][[2]]}

Tableau Formation Code

To use these parameter values to form the Talbeaus we first form the solution by replacing the parameter with the appropriate value. Next, we can automate the process of putting our solution into a Tableau form that will be useable by our odesolvers.  Note that here we have entered a second row of “b’s”. This is so our talbeaus have the correct form to be implemented by ODESolver. The choice of values in the second row only effects the error estimate, so as long as we only used a fixed step driver to test the tableau we should be fine.

solL2=sol/.smallerroptions[[1]];
solLinf=sol/.smallerroptions[[4]];
solrand=sol/.(Subscript[\[FormalC], 2])->RandomReal[{-2,2}];

L2minimize=FormTableau[order,stages,solL2];
Linfminimize=FormTableau[order,stages,solLinf];
rand=FormTableau[order,stages,solrand];

{Map[MatrixForm,L2minimize],Map[MatrixForm,Linfminimize],Map[MatrixForm,rand]}

Results for GF33....

The graph below show the convergence of the three cases. Here the green points are for the parameter chosen for optimality in the L2 norm,, the blue is for the Linf norm and the red is the parameter chosen randomly from the interval [0,2].

a=AsymptoticRKErrorTest[posy0,f,L2minimize,"L2",{1.*(10^(-2)),1.*(10^(-3))},Green];
b=AsymptoticRKErrorTest[posy0,f,Linfminimize,"Linf",{1.*(10^(-2)),1.*(10^(-3))},Blue];
c=AsymptoticRKErrorTest[posy0,f,rand,"rand",{1.*(10^(-2)),1.*(10^(-3))},Red];
Show[a,b,c]

Stability for GF33

Another option for choosing (Subscript[\[FormalC], 2]) would be to choose a value that maximized the stability region? To do this we will form the talbeau maintaining the parameter. The resulting talbeau is:

maxstableTableau=FormTableau[order,stages,sol];
Map[MatrixForm,maxstableTableau]

By doing the standard stability analysis, we generate a function R(z) which is equal to the third order polynomial shown below.

StabilityRegionEquation[maxstableTableau]

The important fact here is that this polynomial is independent of (Subscript[\[FormalC], 2]), thus it doesn’t provide any insight into choosing (Subscript[\[FormalC], 2]). This phenomena is explained in [1, pg 101], and holds whenever the order and stage number of the scheme are the same. Since schemes like that are restriced to having an order four or less, we will look at a larger and more interesting case. Looking at a fifth order six stage scheme, we find that mathematica already has trouble finding an analytic solution using the Solve command. To compensate, we will choose the values of “c” to be evenly spaced values between 0 and 1. Solve then yields a solution of one parameter. The resulting stability function maintains that parameter which implies that different values will yield different stability regions for the resulting 5th order scheme. By trial and error we determine that by setting the parameter (Subscript[\[FormalA], 4,3]) equal to -.165, we get an improved stability region to other comparable schemes. The resulting tableau which we will call GF56 is shown below alongside the accuracy and stability comparision of GF56 (Red) to Fehlberg (blue) and CashKarp (Green).

Tableau creation for GF56

{order,stages}={5,6};
{Eqs,numeqs,numvar}=OrderEqs[order,stages];
var=FormVar[order,stages,numvar];
cons=Table[(Subscript[\[FormalC], i])->(i/stages),{i,2,stages}]
Eqs=Eqs/.cons;

cons=Table[cons[[i]][[1]],{i,1,Length[cons]}]

varnocons=Complement[var,cons]

stabsol=Solve[Eqs,varnocons][[1]]//Quiet

stab=FormTableau[order,stages,stabsol];
Map[MatrixForm,stab]

StabilityRegionEquation[stab]

R[z_,y_]:=1+z+((z^2)/2)+((z^3)/6)+((z^4)/24)+((z^5)/120)-(3/640) (z^6) y

ContourPlot[Norm[R[{x+I y},-.165]]==1,{x,-5,0},{y,-5,5}]

stabtest=stab/.(Subscript[\[FormalA], 4,3])->-.165;

GF56 Results...

Map[MatrixForm,stabtest]

a=AsymptoticRKErrorTest[posy0,f,CashKarp,"stabtest",{1.*(10^(-2)),1.*(10^(-3))},Green];
b=AsymptoticRKErrorTest[posy0,f,Fehlberg,"stabtest",{1.*(10^(-2)),1.*(10^(-3))},Blue];
c=AsymptoticRKErrorTest[posy0,f,stabtest,"stabtest",{1.*(10^(-2)),1.*(10^(-3))},Red];
Show[a,b,c]

a=StabilityRegion[stabtest,Red];
b=StabilityRegion[Fehlberg,Blue];
c=StabilityRegion[CashKarp,Green];
Show[a,b,c]

Higher Order Schemes

As we attempt to create higher order schemes, the number of equations quickly becomes large and difficult to solve. For a 6,8 scheme, even choosing all of the c values still doesn’t allow us to solve the system using the Solve command. To account for this problem we will focus our attention on one specific order equations from each order family[4]. Specifically looking at the equation (∑_i^8)(Subscript[b, i])(Subscript[c, i])=(1/γ), [1, pg 94,95] if we choose the c values we can solve for 6 of the 8 b’s. At this point we can once again use the Solve command to find the nice analytic solution to the system. In this case after making all of our assumptions, our solution is still a function of 4 parameters. To choose these parameters we will once again minimize the principal error terms via the L2 norm. After doing this we arrive at the following scheme called “GF68”. Once created we would like to test the scheme. To do this we will compare its accuracy and stability to Dormand Prince which is the most commonly used 6,8 scheme. The results shown below indicate that our scheme converged slightly faster but did have a smaller stability region.

6th order 8 stage scheme using Solve

Forming the equations

First we set up the equations the same way we did as before.

{order,stages}={6,8};
{Eqs,numeqs,numvar}=OrderEqs[order,stages];
var=FormVar[order,stages,numvar];

lrgmineqs=Table[Eqs[[i]][[1]]-Eqs[[i]][[2]],{i,1,numeqs}];
repl=Table[(Subscript[\[FormalC], i])->(i/stages),{i,2,stages}];
cons=Table[repl[[i]][[1]],{i,1,Length[repl]}];
varnocons=Complement[var,cons];

Solving the b’s and c’s...

repl=Quadsol[order,stages]

Minimizing the Principal Error with four parameters

Eqs=Eqs/.repl;
sol=Solve[Eqs][[1]];

err=Perr[order,stages,sol];
err=err/.repl;
err=err/.sol;

repl=repl/.sol;

sol=Join[sol,repl];

nonfreevar=Join[Table[sol[[i]][[1]],{i,1,Length[sol]}],Table[repl[[i]][[1]],{i,1,Length[repl]}]];
freevar=Complement[varnocons,nonfreevar]

parameters=FindMinimum[Norm[err],freevar]//Quiet

bigsol=sol/.parameters[[2]];

GF68 Results...

GF68=FormTableau[order,stages,bigsol];
GF68=GF68/.parameters[[2]];
Map[MatrixForm,GF68]

a=AsymptoticRKErrorTest[posy0,f,GF68,"GF68",{1.*(10^(-1)),1.*(10^(-2))},Red];
b=AsymptoticRKErrorTest[posy0,f,DormandPrince,"DormandPrince",{1.*(10^(-1)),1.*(10^(-2))},Blue];
Show[{a,b}]

a=StabilityRegion[GF68,Red];
b=StabilityRegion[DormandPrince,Blue];
Show[a,b]

3rd order 3 stage scheme using the minimization technique

Larger/ small Systems

Having run out of options for extending our ability to use the Solve command, we will look at solving the system numerically by re-casting it as a minimization problem. To do this we will subtract the γ(t) [1, pg 97] portion of the order equations to the other side giving us a large nonlinear system of equations all set to be equal zero. We will then minimize (1/2)|ϕ(x)(|^2) . The result might be a solution to ϕ =  0 and thus a solution to our system. As an example, we look at a small 3rd order 3stage problem.

Minimization Equations Code

First we set up the equations the same way we did as before. Next we will create a function using all the information from order equations. We will then perform the minimization as described above. We will do this twice, once with the default tolerances “minsolnonacc” and once with a higher working precision “minsolacc”. To test to see if the minimization does in fact lead to the solution to our system, we will replace the values in the equations with those found in the minimization process. The result should be zero. As shown below, the actual results will depend on the minimization precision.

{order,stages}={3,3};
{Eqs,numeqs,numvar}=OrderEqs[order,stages]

var=FormVar[order,stages,numvar]

mineqs=Table[Eqs[[i]][[1]]-Eqs[[i]][[2]],{i,1,numeqs}]
minsolacc=NMinimize[(1/2) ((Norm[mineqs])^2),var,AccuracyGoal->40,WorkingPrecision->40];
minsolnonacc=NMinimize[(1/2) ((Norm[mineqs])^2),var];

...

mineqs/.minsolacc[[2]]
mineqs/.minsolnonacc[[2]]

Form talbeaus

Next we will form the tableaus that result from the minimization values

minacc=FormTableau[order,stages,minsolacc[[2]]];
minnonacc=FormTableau[order,stages,minsolnonacc[[2]]];

{Map[MatrixForm,minacc],Map[MatrixForm,minnonacc]}

Testing Programs...

After forming the Tableau’s using the minimization process above, we check that the schemes really give us the correct convergence. As we can see in the graph below, both schemes converge like 3rd order, but the scheme formed from the less accurate minimization (green) begins to lose it convergence rate sooner that the scheme formed from the more accurate minimization.

a=AsymptoticRKErrorTest[posy0,f,minacc,"minacc",{1.*(10^(-1)),1.*(10^(-5))},Red];
b=AsymptoticRKErrorTest[posy0,f,minnonacc,"minnonacc",{1.*(10^(-1)),1.*(10^(-5))},Green];
Show[{a,b}]

Large Systems Using minimization techniques

The next step is to take this process and try to find higher order scheme. We attempt to find an 8th order 15 stage scheme. We tried several variations of Nminimize, FindMinimum, and FindInstance, and even some combination such as using a low precision NMinimize to form initial conditions for a FindMinimum. In all cases either the command didn’t finish, we recieved precision errors, or we formed schemes that did not achieve the desired accuracy. In summary we still need either more assumptions or a better technique to handle the large nonlinear systems formed by these order equations.

Setting up the equations

{order,stages}={8,15};
{Eqs,numeqs,numvar}=OrderEqsRow[order,stages];
var=FormVar[order,stages,numvar];

lrgmineqs=Table[Eqs[[i]][[1]]-Eqs[[i]][[2]],{i,1,numeqs}];
repl=Table[(Subscript[\[FormalC], i])->(i/stages),{i,2,stages}];
cons=Table[repl[[i]][[1]],{i,1,Length[repl]}];
varnocons=Complement[var,cons];

repl=Quadsol[order,stages];

lrgmineqs=lrgmineqs/.repl;

nonfreevar=Table[repl[[i]][[1]],{i,1,Length[repl]}];
freevar=Complement[varnocons,nonfreevar];

Results

Timing[lrgminsolacc=NMinimize[(1/2) ((Norm[lrgmineqs])^2),freevar,(*AccuracyGoal->20,WorkingPrecision->30,*)MaxIterations->100];]

init=Table[{lrgminsolacc[[2]][[i]][[1]],lrgminsolacc[[2]][[i]][[2]]},{i,1,Length[lrgminsolacc[[2]]]}];

Timing[findminsol=FindMinimum[(1/2) ((Norm[lrgmineqs])^2),init,AccuracyGoal->30,WorkingPrecision->60,MaxIterations->100];]

lrgminsolacc[[1]]

....

lrgmineqs/.lrgminsolacc[[2]]//N

lrgmintab=FormTableau[order,stages,lrgminsolacc[[2]]];

References

[1] Butcher J. C. (2008). Numerical Methods for Ordinary Differential Equations. Wiley
[2] Butcher J.C. (). The non-existence of ten stage eighth order explicit Runge-Kutta methods. BIT Numerical Mathematics, 25 (3), 521-540
[3] http://reference.wolfram.com/mathematica/NumericalDifferentialEquationAnalysis/tutorial/NumericalDifferentialEquationAnalysis.html 
[4] Butcher J.C. (2008). http://www.math.auckland.ac.nz/~butcher/ODE-book-2008/Tutorials/IRK.pdf