update tests; closes #44

README.md

@@ -23,6 +23,9 @@ A common workflow for designing Runge-Kutta methods is to use **polyopt** to fin
 appropriate stability function and then **RK-coeff-opt** to determine the Runge-Kutta
 method coefficients.
 
+To run the tests, execute the MATLAB script `test.m`.
+
+
 # Authors
 The code is primarily developed and maintained by David Ketcheson.
 The following people have also made important constributions to RK-opt (listed alphabetically):

RK-coeff-opt/README.rst

@@ -1,4 +1,10 @@
-This subpackage contains routines for finding optimal Runge-Kutta method coefficents,
-given a prescribed order of accuracy, number of stages, and an objective function.  
+This subpackage contains routines for finding optimal Runge-Kutta method coefficients,
+given a prescribed order of accuracy, number of stages, and an objective function.
 Constraints on the stability polynomial (possibly obtained using **polyopt** or **am_radius-opt**)
 can optionally be provided.
+
+To run the tests, execute the MATLAB commands
+```
+results_rkopt = runtests('test_rkopt.m');
+table(results_rkopt)
+```

RK-coeff-opt/oc_ksrk.m

@@ -2,7 +2,7 @@
 % function coneq= oc_ksrk(A,b,D,theta,p)
 % Order conditions for multistep-RK methods.
 %
-% ..warning:: 
+% ..warning::
 %
 %         Here we assume a certain minimum stage order,
 %         which is necessarily true for methods with
@@ -40,21 +40,24 @@
 
 if p>=2
     coneq(2)=b'*c   -(1-(-ll').^2*theta')/2;
+end
 if p>=3
     coneq(3)=b'*c.^2-(1-(-ll').^3*theta')/3;
     coneq(4)=b'*tau_2;
+end
 if p>=4
     coneq(5)=b'*c.^3-(1-(-ll').^4*theta')/4;
     coneq(6)=b'*C*tau_2;
     coneq(7)=b'*A*tau_2;
     coneq(8)=b'*tau_3;
+end
 if p==5          % for p>=5 we assume tau_2=0
     coneq(9)=b'*c.^4-(1-(-ll').^5*theta')/5;
     coneq(10)=b'*A*tau_3;
     coneq(11)=b'*tau_4;
     coneq(12)=b'*C*tau_3;
     coneq=[coneq tau_2'];
-elseif p==6 
+elseif p==6
     coneq(9)=b'*c.^4-(1-(-ll').^5*theta')/5;
     coneq(10)=b'*A*tau_3;
     coneq(11)=b'*tau_4;
@@ -166,8 +169,8 @@
     coneq(5)  = b'*c.^4-(1-(-ll').^5*theta')/5;
     coneq(6)  = b'*c.^5-(1-(-ll').^6*theta')/6;
     coneq(7)  = b'*c.^6-(1-(-ll').^7*theta')/7;
-    coneq(8)  = b'*c.^7-(1-(-ll').^8*theta')/8;	
-    coneq(9)  = b'*c.^8-(1-(-ll').^9*theta')/9;    
+    coneq(8)  = b'*c.^7-(1-(-ll').^8*theta')/8;
+    coneq(9)  = b'*c.^8-(1-(-ll').^9*theta')/9;
     coneq(10) = b'*A*tau_5;
     coneq(11) = b'*tau_6;
     coneq(12) = b'*C*tau_5;
@@ -195,7 +198,7 @@
     coneq(34) = b'*A*C*tau_6;
     coneq(35) = b'*A*C^2*tau_5;
     coneq(36) = b'*tau_8;
-    coneq(37) = b'*c.^9-(1-(-ll').^10*theta')/10; 
+    coneq(37) = b'*c.^9-(1-(-ll').^10*theta')/10;
     coneq(38) = b'*A^2*tau_5;
     coneq(39) = b'*A*tau_6;
     coneq(40) = b'*A*C*tau_5;
@@ -253,29 +256,29 @@
     coneq(92) = b'*tau_9;
 
     coneq=[coneq tau_2' tau_3' tau_4'];
-  
+
 elseif p==11; % Assume stage order 5
     coneq(2)  = b'*c-(1-(-ll').^2*theta')/2;
     coneq(3)  = b'*c.^2-(1-(-ll').^3*theta')/3;
     coneq(4)  = b'*c.^3-(1-(-ll').^4*theta')/4;
     coneq(5)  = b'*c.^4-(1-(-ll').^5*theta')/5;
     coneq(6)  = b'*c.^5-(1-(-ll').^6*theta')/6;
     coneq(7)  = b'*c.^6-(1-(-ll').^7*theta')/7;
-    coneq(8)  = b'*c.^7-(1-(-ll').^8*theta')/8;	
-    coneq(9)  = b'*c.^8-(1-(-ll').^9*theta')/9;      
+    coneq(8)  = b'*c.^7-(1-(-ll').^8*theta')/8;
+    coneq(9)  = b'*c.^8-(1-(-ll').^9*theta')/9;
     coneq(10) = b'*c.^9-(1-(-ll').^10*theta')/10;
     coneq(11) = b'*c.^10-(1-(-ll').^11*theta')/11;
 
     coneq(12) = b'*tau_6;
     coneq(13) = b'*A*tau_6;
-    coneq(14) = b'*tau_7;     
+    coneq(14) = b'*tau_7;
     coneq(15) = b'*C*tau_6;
     coneq(16) = b'*C*A*tau_6;
     coneq(17) = b'*C*tau_7;
     coneq(18) = b'*C^2*tau_6;
     coneq(19) = b'*A^2*tau_6;
-    coneq(20) = b'*A*tau_7;       
-    coneq(21) = b'*A*C*tau_6;       
+    coneq(20) = b'*A*tau_7;
+    coneq(21) = b'*A*C*tau_6;
     coneq(22) = b'*tau_8;
     coneq(23) = b'*A*tau_6;
     coneq(24) = b'*A^2*tau_6;
@@ -303,14 +306,14 @@
 
     coneq(46) = b'*A*tau_6;
     coneq(47) = b'*A*A*tau_6;
-    coneq(48) = b'*A*tau_7;     
+    coneq(48) = b'*A*tau_7;
     coneq(49) = b'*A*C*tau_6;
     coneq(50) = b'*A*C*A*tau_6;
     coneq(51) = b'*A*C*tau_7;
     coneq(52) = b'*A*C^2*tau_6;
     coneq(53) = b'*A*A^2*tau_6;
-    coneq(54) = b'*A*A*tau_7;       
-    coneq(55) = b'*A*A*C*tau_6;       
+    coneq(54) = b'*A*A*tau_7;
+    coneq(55) = b'*A*A*C*tau_6;
     coneq(56) = b'*A*tau_8;
     coneq(57) = b'*A*A*tau_6;
     coneq(58) = b'*A*A^2*tau_6;
@@ -334,18 +337,18 @@
     coneq(76) = b'*A*C*A*tau_7;
     coneq(77) = b'*A*C*A*C*tau_6;
     coneq(78) = b'*A*C*tau_8;
-    coneq(79) = b'*A*tau_9;     
+    coneq(79) = b'*A*tau_9;
 
     coneq(80) = b'*C*tau_6;
     coneq(81) = b'*C*A*tau_6;
-    coneq(82) = b'*C*tau_7;     
+    coneq(82) = b'*C*tau_7;
     coneq(83) = b'*C*C*tau_6;
     coneq(84) = b'*C*C*A*tau_6;
     coneq(85) = b'*C*C*tau_7;
     coneq(86) = b'*C*C^2*tau_6;
     coneq(87) = b'*C*A^2*tau_6;
-    coneq(88) = b'*C*A*tau_7;       
-    coneq(89) = b'*C*A*C*tau_6;       
+    coneq(88) = b'*C*A*tau_7;
+    coneq(89) = b'*C*A*C*tau_6;
     coneq(90) = b'*C*tau_8;
     coneq(91) = b'*C*A*tau_6;
     coneq(92) = b'*C*A^2*tau_6;
@@ -371,9 +374,9 @@
     coneq(112) = b'*C*C*tau_8;
     coneq(113) = b'*C*tau_9;
     coneq(114) = b'*tau_10;
-        
+
     coneq=[coneq tau_2' tau_3' tau_4' tau_5'];
 
-elseif p>11 
+elseif p>11
     disp('Order conditions for p>11 are not coded up yet');
 end