(*-------------------- algorithm for derivation of Runge-Kutta pairs of orders 8(7) --------------------*) c2=1/16;c5=39/100;c6=7/15;c7=39/250;c8=24/25;c10=11/12;c11=19/20;a87=1993321838240/380523459069; ww12=0;ww13=-45884771325/82380821152;w13=0; (*DVERK78*) c9=(14*c6^2*(7*c7^2*c8+c7*(7*c8^2-12*c8+1)+c8)+c6*(14*c7^2*(7*c8^2-12*c8+1)-7*c7*(24*c8^2-33*c8+4)+14*c8^2 -28*c8+3)+14*c7^2*c8+c7*(14*c8^2-28*c8+3)+3*c8)/(2*(7*c6^2*(7*c7^2*(15*c8^2-10*c8+2)-2*c7*(35*c8^2 -26*c8+6)+14*c8^2-12*c8+3)-7*c6*(2*c7^2*(35*c8^2-26*c8+6)-c7*(52*c8^2-42*c8+11)+12*c8^2-11*c8+3) +7*c7^2*(14*c8^2-12*c8+3)-7*c7*(12*c8^2-11*c8+3)+21*c8^2-21*c8+6)); w12=-(840*w13*(c9-1)*(c8-1)*(c7-1)*(c6-1)*(c11-1)*(c10-1)-2*c10*(7*c11*(c6*(5*c7*(2*c8*(3*c9-2)-4*c9+3) -5*c8*(4*c9-3)+3*(5*c9-4))+c7*(3*(5*c9-4)-5*c8*(4*c9-3))+3*c8*(5*c9-4)-2*(6*c9-5)) -7*c6*(c7*(5*c8*(4*c9-3)-3*(5*c9-4))-3*c8*(5*c9-4)+2*(6*c9-5))+7*c7*(3*c8*(5*c9-4)-2*(6*c9-5)) -2*(7*c8*(6*c9-5)-5*(7*c9-6)))+2*c11*(7*c6*(c7*(5*c8*(4*c9-3)-3*(5*c9-4))-3*c8*(5*c9-4)+2*(6*c9-5)) -7*c7*(3*c8*(5*c9-4)-2*(6*c9-5))+2*(7*c8*(6*c9-5)-5*(7*c9-6)))-2*c6*(7*c7*(3*c8*(5*c9-4)-2*(6*c9-5)) -2*(7*c8*(6*c9-5)-5*(7*c9-6)))+4*c7*(7*c8*(6*c9-5)-5*(7*c9-6))-5*(4*c8*(7*c9-6)-3*(8*c9-7)))/ (840*(c10-1)*(c11-1)*(c6-1)*(c7-1)*(c8-1)*(c9-1)); w11=-(420*w12*(c9-1)*(c8-1)*(c7-1)*(c6-1)*(c10-1)+420*w13*(c9-1)*(c8-1)*(c7-1)*(c6-1)*(c10-1) -7*c10*(c6*(5*c7*(2*c8*(3*c9-2)-4*c9+3)-5*c8*(4*c9-3)+3*(5*c9-4))+c7*(3*(5*c9-4)-5*c8*(4*c9-3)) +3*c8*(5*c9-4)-2*(6*c9-5))+7*c6*(c7*(5*c8*(4*c9-3)-3*(5*c9-4))-3*c8*(5*c9-4)+2*(6*c9-5)) -7*c7*(3*c8*(5*c9-4)-2*(6*c9-5))+2*(7*c8*(6*c9-5)-5*(7*c9-6)))/ (420*c11*(c10-c11)*(c11-c6)*(c11-c7)*(c11-c8)*(c11-c9)); w10=-(60*w11*c11*(c11^4-c11^3*(c6+c7+c8+c9)+c11^2*(c6*(c7+c8+c9)+c7*(c8+c9)+c8*c9)-c11*(c6*(c7*(c8+c9) +c8*c9)+c7*c8*c9)+c6*c7*c8*c9)+60*w12*(c9-1)*(c8-1)*(c7-1)*(c6-1)+60*w13*(c9-1)*(c8-1)*(c7-1)*(c6-1) -c6*(5*c7*(2*c8*(3*c9-2)-4*c9+3)-5*c8*(4*c9-3)+3*(5*c9-4))+c7*(5*c8*(4*c9-3)-3*(5*c9-4)) -3*c8*(5*c9-4)+2*(6*c9-5))/(60*c10*(c10-c6)*(c10-c7)*(c10-c8)*(c10-c9)); w9=(60*w10*c10*(c10^3-c10^2*(c6+c7+c8)+c10*(c6*(c7+c8)+c7*c8)-c6*c7*c8)+60*w11*c11*(c11^3-c11^2*(c6+c7+c8) +c11*(c6*(c7+c8)+c7*c8)-c6*c7*c8)-60*w12*(c8-1)*(c7-1)*(c6-1)-60*w13*(c8-1)*(c7-1)*(c6-1) +5*c6*(2*c7*(3*c8-2)-4*c8+3)-5*c7*(4*c8-3)+3*(5*c8-4))/(60*c9*(c6-c9)*(c7-c9)*(c8-c9)); w8=-(12*w10*c10*(c10^2-c10*(c6+c7)+c6*c7)+12*w11*c11*(c11^2-c11*(c6+c7)+c6*c7)+12*w12*(c7-1)*(c6-1) +12*w13*(c7-1)*(c6-1)+12*w9*c9*(c7-c9)*(c6-c9)-2*c6*(3*c7-2)+4*c7-3)/(12*c8*(c6-c8)*(c7-c8)); w7=(6*w10*c10*(c10-c6)+6*w11*c11*(c11-c6)-6*w12*(c6-1)-6*w13*(c6-1)-6*w8*c8*(c6-c8)-6*w9*c9*(c6-c9)+3*c6-2)/ (6*c7*(c6-c7)); w6=-(2*w10*c10+2*w11*c11+2*w12+2*w13+2*w7*c7+2*w8*c8+2*w9*c9-1)/(2*c6); ww11=-(420*ww12*(c9-1)*(c8-1)*(c7-1)*(c6-1)*(c10-1)+420*ww13*(c9-1)*(c8-1)*(c7-1)*(c6-1)*(c10-1) -7*c10*(c6*(5*c7*(2*c8*(3*c9-2)-4*c9+3)-5*c8*(4*c9-3)+3*(5*c9-4))+c7*(3*(5*c9-4)-5*c8*(4*c9-3))+3*c8*(5*c9-4) -2*(6*c9-5))+7*c6*(c7*(5*c8*(4*c9-3)-3*(5*c9-4))-3*c8*(5*c9-4)+2*(6*c9-5))-7*c7*(3*c8*(5*c9-4)-2*(6*c9-5)) +2*(7*c8*(6*c9-5)-5*(7*c9-6)))/(420*c11*(c10-c11)*(c11-c6)*(c11-c7)*(c11-c8)*(c11-c9)); ww10=-(60*ww11*c11*(c11^4-c11^3*(c6+c7+c8+c9)+c11^2*(c6*(c7+c8+c9)+c7*(c8+c9)+c8*c9)-c11*(c6*(c7*(c8+c9)+c8*c9)+c7*c8*c9) +c6*c7*c8*c9)+60*ww12*(c9-1)*(c8-1)*(c7-1)*(c6-1)+60*ww13*(c9-1)*(c8-1)*(c7-1)*(c6-1)-c6*(5*c7*(2*c8*(3*c9-2)-4*c9+3) -5*c8*(4*c9-3)+3*(5*c9-4))+c7*(5*c8*(4*c9-3)-3*(5*c9-4))-3*c8*(5*c9-4)+2*(6*c9-5))/ (60*c10*(c10-c6)*(c10-c7)*(c10-c8)*(c10-c9)); ww9=(60*ww10*c10*(c10^3-c10^2*(c6+c7+c8)+c10*(c6*(c7+c8)+c7*c8)-c6*c7*c8)+60*ww11*c11*(c11^3-c11^2*(c6+c7+c8) +c11*(c6*(c7+c8)+c7*c8)-c6*c7*c8)-60*ww12*(c8-1)*(c7-1)*(c6-1)-60*ww13*(c8-1)*(c7-1)*(c6-1)+5*c6*(2*c7*(3*c8-2)-4*c8+3) -5*c7*(4*c8-3)+3*(5*c8-4))/(60*c9*(c6-c9)*(c7-c9)*(c8-c9)); ww8=-(12*ww10*c10*(c10^2-c10*(c6+c7)+c6*c7)+12*ww11*c11*(c11^2-c11*(c6+c7)+c6*c7)+12*ww12*(c7-1)*(c6-1)+12*ww13*(c7-1)*(c6-1) +12*ww9*c9*(c7-c9)*(c6-c9)-2*c6*(3*c7-2)+4*c7-3)/(12*c8*(c6-c8)*(c7-c8)); ww7=(6*ww10*c10*(c10-c6)+6*ww11*c11*(c11-c6)-6*ww12*(c6-1)-6*ww13*(c6-1)-6*ww8*c8*(c6-c8)-6*ww9*c9*(c6-c9)+3*c6-2)/(6*c7*(c6-c7)); ww6=-(2*ww10*c10+2*ww11*c11+2*ww12+2*ww13+2*ww7*c7+2*ww8*c8+2*ww9*c9-1)/(2*c6); a1110=-(2*c6*(7*c7*(c8*(5*c9-2)-2*c9+1)-7*c8*(2*c9-1)+7*c9-4)-2*c7*(7*c8*(2*c9-1)-7*c9+4)+2*c8*(7*c9-4)-8*c9+5)/ (1680*w11*c10*(c10-c6)*(c10-c7)*(c10-c8)*(c10-c9)*(c11-1)); a109=(3*c11*(7*c6*(c7*(5*c8-2)-2*c8+1)-7*c7*(2*c8-1)+7*c8-4)-c6*(7*c7*(9*c8-4)-28*c8+15)+c7*(28*c8-15)-3*(5*c8-3))/ (2520*w10*c9*(c10-1)*(c10-c11)*(c6-c9)*(c7-c9)*(c8-c9)); a119=-(1680*a109*w10*c9*(c8-c9)*(c7-c9)*(c6-c9)*(c10^2-c10*(c9+1)+c9)+2*c10*(7*c6*(c7*(5*c8-2)-2*c8+1) -7*c7*(2*c8-1)+7*c8-4)-2*c6*(7*c7*(2*c8-1)-7*c8+4)+2*c7*(7*c8-4)-8*c8+5)/ (1680*w11*c9*(c10-c9)*(c11-1)*(c6-c9)*(c7-c9)*(c8-c9)); c4=c6*(4*c5-3*c6)/(2*(3*c5-2*c6)); c3=2*c4/3; a32=c3^2/(2*c2); a43=c4^2/(2*c3); a54=c5^2*(3*c3-2*c5)/(6*c4*(c3-c4)); a53=-(2*a54*c4-c5^2)/(2*c3); a65=c6^2*(3*c4-2*c6)/(6*c5*(c4-c5)); a64=-(2*a65*c5-c6^2)/(2*c4); a76=c7^2*(2*c4*(3*c5-2*c7)-c7*(4*c5-3*c7))/(12*c6*(c4-c6)*(c5-c6)); a75=-(6*a76*c6*(c4-c6)-c7^2*(3*c4-2*c7))/(6*c5*(c4-c5)); a74=-(2*a75*c5+2*a76*c6-c7^2)/(2*c4); a86=-(12*a87*c7*(c5-c7)*(c4-c7)-c8^2*(2*c4*(3*c5-2*c8)-c8*(4*c5-3*c8)))/(12*c6*(c4-c6)*(c5-c6)); a85=-(6*a86*c6*(c4-c6)+6*a87*c7*(c4-c7)-c8^2*(3*c4-2*c8))/(6*c5*(c4-c5)); a84=-(2*a85*c5+2*a86*c6+2*a87*c7-c8^2)/(2*c4); a114=-(a64*(c6-1)*(w10*(w6*ww9*(c10-c6)+w9*ww6*(c9-c10))+w6*w9*ww10*(c6-c9))+a74*(c7-1)*(w10*(w7*ww9*(c10-c7) +w9*ww7*(c9-c10))+w7*w9*ww10*(c7-c9))+a84*(c8-1)*(w10*(w8*ww9*(c10-c8)+w9*ww8*(c9-c10))+w8*w9*ww10*(c8-c9)))/ ((w10*(w11*ww9*(c10-c11)+w9*ww11*(c9-c10))+w11*w9*ww10*(c11-c9))*(c11-1)); a104=-(a114*w11*(c11^2-c11*(c9+1)+c9)+a64*w6*(c6^2-c6*(c9+1)+c9)+a74*w7*(c7^2-c7*(c9+1)+c9)+a84*w8*(c8^2-c8*(c9+1)+c9))/ (w10*(c10-1)*(c10-c9)); a94=(a104*w10*(c10^2-c10*(c11+1)+c11)+a64*w6*(1-c6)*(c11-c6)+a74*w7*(1-c7)*(c11-c7)+a84*w8*(1-c8)*(c11-c8))/(w9*(c11-c9)*(c9-1)); a115=-(a65*(c6-1)*(w10*(w6*ww9*(c10-c6)+w9*ww6*(c9-c10))+w6*w9*ww10*(c6-c9))+a75*(c7-1)*(w10*(w7*ww9*(c10-c7) +w9*ww7*(c9-c10))+w7*w9*ww10*(c7-c9))+a85*(c8-1)*(w10*(w8*ww9*(c10-c8)+w9*ww8*(c9-c10))+w8*w9*ww10*(c8-c9)))/ ((w10*(w11*ww9*(c10-c11)+w9*ww11*(c9-c10))+w11*w9*ww10*(c11-c9))*(c11-1)); a105=-(a115*w11*(c11^2-c11*(c9+1)+c9)+a65*w6*(c6^2-c6*(c9+1)+c9)+a75*w7*(c7^2-c7*(c9+1)+c9)+a85*w8*(c8^2-c8*(c9+1)+c9))/ (w10*(c10-1)*(c10-c9)); a95=(a105*w10*(c10^2-c10*(c11+1)+c11)+a65*w6*(1-c6)*(c11-c6)+a75*w7*(1-c7)*(c11-c7)+a85*w8*(1-c8)*(c11-c8))/(w9*(c11-c9)*(c9-1)); a98=-(12*a94*c4*(c4^2-c4*(c6+c7)+c6*c7)+12*a95*c5*(c5^2-c5*(c6+c7)+c6*c7)-c9^2*(2*c6*(3*c7-2*c9)-c9*(4*c7-3*c9)))/ (12*c8*(c6-c8)*(c7-c8)); a97=(6*a94*c4*(c4-c6)+6*a95*c5*(c5-c6)-6*a98*c8*(c6-c8)+c9^2*(3*c6-2*c9))/(6*c7*(c6-c7)); a96=-(2*a94*c4+2*a95*c5+2*a97*c7+2*a98*c8-c9^2)/(2*c6); a108=-(12*a104*c4*(c4^2-c4*(c6+c7)+c6*c7)+12*a105*c5*(c5^2-c5*(c6+c7)+c6*c7)+12*a109*c9*(c7-c9)*(c6-c9)-c10^2*(3*c10^2 -4*c10*(c6+c7)+6*c6*c7))/(12*c8*(c6-c8)*(c7-c8)); a107=(6*a104*c4*(c4-c6)+6*a105*c5*(c5-c6)-6*a108*c8*(c6-c8)-6*a109*c9*(c6-c9)-c10^2*(2*c10-3*c6))/(6*c7*(c6-c7)); a106=-(2*a104*c4+2*a105*c5+2*a107*c7+2*a108*c8+2*a109*c9-c10^2)/(2*c6); a118=-(12*a1110*c10*(c10^2-c10*(c6+c7)+c6*c7)+12*a114*c4*(c4^2-c4*(c6+c7)+c6*c7)+12*a115*c5*(c5^2-c5*(c6+c7)+c6*c7) +12*a119*c9*(c7-c9)*(c6-c9)-c11^2*(3*c11^2-4*c11*(c6+c7)+6*c6*c7))/(12*c8*(c6-c8)*(c7-c8)); a117=(6*a1110*c10*(c10-c6)+6*a114*c4*(c4-c6)+6*a115*c5*(c5-c6)-6*a118*c8*(c6-c8)-6*a119*c9*(c6-c9)-c11^2*(2*c11-3*c6))/(6*c7*(c6-c7)); a116=-(2*a1110*c10+2*a114*c4+2*a115*c5+2*a117*c7+2*a118*c8+2*a119*c9-c11^2)/(2*c6); a1311=(c11-1)*(w11*ww12-w12*ww11)/(w12*ww13-w13*ww12); a1211=-(a1311*w13+w11*(c11-1))/w12; a1310=(a1110*(w11*ww12-w12*ww11)+(c10-1)*(w10*ww12-w12*ww10))/(w12*ww13-w13*ww12); a1210=-(a1110*w11+a1310*w13+w10*(c10-1))/w12; a139=(a109*(w10*ww12-w12*ww10)+a119*(w11*ww12-w12*ww11)+(1-c9)*(w12*ww9-w9*ww12))/(w12*ww13-w13*ww12); a129=-(a109*w10+a119*w11+a139*w13+w9*(c9-1))/w12; a138=(a108*(w10*ww12-w12*ww10)+a118*(w11*ww12-w12*ww11)+a98*(w9*ww12-w12*ww9)+(1-c8)*(w12*ww8-w8*ww12))/(w12*ww13-w13*ww12); a128=-(a108*w10+a118*w11+a138*w13+a98*w9+w8*(c8-1))/w12; a137=(a107*(w10*ww12-w12*ww10)+a117*(w11*ww12-w12*ww11)+a87*(w8*ww12-w12*ww8)+a97*(w9*ww12-w12*ww9)+(1-c7)*(w12*ww7-w7*ww12))/ (w12*ww13-w13*ww12); a127=-(a107*w10+a117*w11+a137*w13+a87*w8+a97*w9+w7*(c7-1))/w12; a136=(a106*(w10*ww12-w12*ww10)+a116*(w11*ww12-w12*ww11)+a76*(w7*ww12-w12*ww7)+a86*(w8*ww12-w12*ww8)+a96*(w9*ww12-w12*ww9) +(1-c6)*(w12*ww6-w6*ww12))/(w12*ww13-w13*ww12); a126=-(a106*w10+a116*w11+a136*w13+a76*w7+a86*w8+a96*w9+w6*(c6-1))/w12; a135=(a105*(w10*ww12-w12*ww10)+a115*(w11*ww12-w12*ww11)+a65*(w6*ww12-w12*ww6)+a75*(w7*ww12-w12*ww7)+a85*(w8*ww12-w12*ww8) +a95*(w9*ww12-w12*ww9))/(w12*ww13-w13*ww12); a125=-(a105*w10+a115*w11+a135*w13+a65*w6+a75*w7+a85*w8+a95*w9)/w12; a134=(a104*(w10*ww12-w12*ww10)+a114*(w11*ww12-w12*ww11)+a64*(w6*ww12-w12*ww6)+a74*(w7*ww12-w12*ww7) +a84*(w8*ww12-w12*ww8)+a94*(w9*ww12-w12*ww9))/(w12*ww13-w13*ww12); a124=-(a104*w10+a114*w11+a134*w13+a64*w6+a74*w7+a84*w8+a94*w9)/w12; a21=c2; a31=c3-a32; a41=c4-a43; a51=-a53-a54+c5; a61=-a64-a65+c6; a71=-a74-a75-a76+c7; a81=-a84-a85-a86-a87+c8; a91=-a94-a95-a96-a97-a98+c9; a101=-a104-a105-a106-a107-a108-a109+c10; a111=-a1110-a114-a115-a116-a117-a118-a119+c11; a121=-a1210-a1211-a124-a125-a126-a127-a128-a129+1; a131=-a1310-a1311-a134-a135-a136-a137-a138-a139+1; c={0,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,1,1}; w={1-w6-w7-w8-w9-w10-w11-w12-w13,0,0,0,0,w6,w7,w8,w9,w10,w11,w12,w13}; ww={1-ww6-ww7-ww8-ww9-ww10-ww11-ww12-ww13,0,0,0,0,ww6,ww7,ww8,ww9,ww10,ww11,ww12,ww13}; e={1,1,1,1,1,1,1,1,1,1,1,1,1}; a={{0,0,0,0,0,0,0,0,0,0,0,0,0}, {a21,0,0,0,0,0,0,0,0,0,0,0,0}, {a31,a32,0,0,0,0,0,0,0,0,0,0,0}, {a41,0,a43,0,0,0,0,0,0,0,0,0,0}, {a51,0,a53,a54,0,0,0,0,0,0,0,0,0}, {a61,0,0,a64,a65,0,0,0,0,0,0,0,0}, {a71,0,0,a74,a75,a76,0,0,0,0,0,0,0}, {a81,0,0,a84,a85,a86,a87,0,0,0,0,0,0}, {a91,0,0,a94,a95,a96,a97,a98,0,0,0,0,0}, {a101,0,0,a104,a105,a106,a107,a108,a109,0,0,0,0}, {a111,0,0,a114,a115,a116,a117,a118,a119,a1110,0,0,0}, {a121,0,0,a124,a125,a126,a127,a128,a129,a1210,a1211,0,0}, {a131,0,0,a134,a135,a136,a137,a138,a139,a1310,a1311,0,0}}; (*-------------------------------------------------------------------------------------------------------------------------------*) (*------------------------------- Runge-Kutta pair of orders 8(7) with extented stability region --------------------------------*) a={{0,0,0,0,0,0,0,0,0,0,0,0,0}, {0.0100000000000548001,0,0,0,0,0,0,0,0,0,0,0,0}, {-0.361866893966954953,0.457525156145154926,0,0,0,0,0,0,0,0,0,0,0}, {0.0358718483168249898,0,0.107615544950474969,0,0,0,0,0,0,0,0,0,0}, {0.491468336121696283,0,-1.88065074984986715,1.77649227042813523,0,0,0,0,0,0,0,0,0}, {0.0390372212808226731,0,0,0.230729570202058204,0.202724163763739067,0,0,0,0,0,0,0,0}, {0.0566763049610813938,0,0,0.113539003895125816,-0.0239102037832446040,0.0114581245323579544,0,0,0,0,0,0,0}, {1.38561312043069786,0,0,-62.7210606385233894,-13.7145038445302680,8.84040388475592411,67.1716035927809543,0,0,0,0,0,0}, {0.882621875282926857,0,0,-40.3854408272919848,-9.25083243548891832,6.21178463330449754,43.4502111365405483, 0.00865193799487277629,0,0,0,0,0}, {0.476973554680081405,0,0,-21.5445578307502176,-5.22838516660087428,3.76899052206361569,23.3654344382504739, -0.0274025414042129872,0.0614792745518951220,0,0,0,0}, {2.10260074710297510,0,0,-92.8363996221548944,-19.1136725662967247,11.8143320479647333,98.9910922797547281,0.173025655670426156, -0.571919557178592258,0.457863161333586103,0,0,0}, {1.10869611538397427,0,0,-47.9518216827019317,-9.59835412423178996,5.95042502117200074,51.2234610576634570,0.848435964263701314, -2.42300586617722244,1.81712170785586112,0.0250418067719496143,0,0}, {2.46427522909090341,0,0,-110.847673536399517,-23.6425789507075051,14.8559372997165122,118.328016591970549,-0.594981312329926486, 1.50683777297286678,-1.06461201636364725,-0.00522107795023598360,0,0}}; w={0.0450940925252529857,0,0,0,0,0.358463860936017270,0.252236408588385151,-0.676822967127669664,-0.787174860970923295, 0.999999985482866741,-1.00000057891919183,0.871130778817386403,0.937073280667876240}; ww={0.04507973308293407,0,0,0,0,0.35823849913486688,0.25228983678114392,-0.41344736135123274,-0.95009726718876393, 1.04583916863402738,-0.86850949386135742,0.74971248421181832,0.78089440055656353}; c={0,0.0100000000000548001,0.0956582621781999727,0.143487393267299959,0.387309856699964361,0.472490955246619944, 0.157763229605320560,0.962056114913918869,0.916996320341942389,0.872532250790761222,1.01692214619623742,1,1}; (*-------------------------------------------------------------------------------------------------------------------------------*)