# # 2019-08-30 List of Runge-Kutta methods # (c) 2002-2019 Sergey Khashin # e-mail: khash2@gmail.com # www: http://math.ivanovo.ac.ru/dalgebra/Khashin/index.html # # ---- RK-methods --------------------------------------- # RK_44 (c2, c3) # RK_44a(b3) # RK_44b(b4) # RK_45 (c2, c3, c4, c5, b2, a32, a42) # RK_56 (c2, c3, c4, c5, a43) b2 = 0 # RK_56a(c2, c3, c4, c5, b2, a43) b2 != 0 # RK_67 (c2, c3, c5, c6) b2 =0 # RK_67a(c3, c5, c6) b2<>0 # RK_67b(c2, a32, b3) c2=c3=c6 # RK_79(c4,c5) from Butchers book # RK_79(c4,c5) from Butchers book # RK_8_11(a105,b8) from Assui Kouassi Richard, Seka Hippolyte # --------------------------------------------------------------------- # # # RK(4,4), example: # a4 := RK_44(1/3, 1/2); RK_44 := proc(c2, c3)::Matrix; local a,d; if c2=1/2 then ERROR(`RK_44: c2=1/2`); end if; if c2=c3 then ERROR(`RK_44: c2=c3`); end if; if c2=1 then ERROR(`RK_44: c2=1`); end if; if c3=1 then ERROR(`RK_44: c3=1`); end if; d := factor(3-4*(c2+c3)+6*c2*c3); if d=0 then ERROR(`RK_set44: 3-4*(c2+c3)+6*c2*c3=0`); end if; a := Matrix(5,5); a[2,1] := c2; a[3,2] := c3*(c3-c2)/(2*c2*(1-2*c2)); a[3,1] := c3-a[3,2]; a[4,2] := (c2-1)*(2-c2-5*c3+4*c3*c3)/(2*c2*(c3-c2)*d); a[4,3] := (1-2*c2)*(c3-1)*(c2-1)/(c3*(c3-c2)*d); a[4,1] := 1 - a[4,3] - a[4,2]; a[5,2] := (1-2*c3)/(12*c2*(c2-1)*(c3-c2)); a[5,3] := (2*c2-1)/(12*c3*(c3-1)*(c3-c2)); a[5,4] := (3-4*c2-4*c3+6*c2*c3)/(12*(c2-1)*(c3-1)); a[5,1] := 1 - a[5,4] - a[5,3] - a[5,2]; simplify(a); end proc; # --------------------------------------------------------------------- # RK(4,4) c2=c3=1/2: RK_44a := proc(b3)::Matrix; local a; if b3=0 then ERROR(`RK_44a: b3=0`); end if; a := Matrix(5,5); a[2,1] := 1/2; a[3,2] := 1/6/b3; a[3,1] := 1/2-a[3,2]; a[4,2] := 1-3*b3; a[4,3] := 3*b3; a[4,1] := 0; a[5,2] := 2/3-b3; a[5,3] := b3; a[5,4] := 1/6; a[5,1] := 1/6; simplify(a); end proc; # --------------------------------------------------------------------- # RK(4,4) c2=1, c3=1/2: RK_44b := proc(b4)::Matrix; local a; if b4=0 then ERROR(`RK_44b: b4=0`); end if; a := Matrix(5,5); a[2,1] := 1; a[3,1] := 3/8; a[3,2] := 1/8; a[4,3] := 1/3/b4; a[4,2] := -a[4,3]/4; a[4,1] := 1-a[4,2]-a[4,3]; a[5,1] := 1/6; a[5,2] := 1/6-b4; a[5,3] := 2/3; a[5,4] := b4; simplify(a); end proc; # --------------------------------------------------------------------- RK_45 := proc(c2, c3, c4, c5, b2, a32, a42)::Matrix; local a,d,d1; a := Matrix(6,6); if c3=0 then ERROR(`RK_45: c3=0`); end if; d := factor(c3^2-c5*c3+c5*c4-c4*c3); if d =0 then ERROR(`RK_45: c3^2-c5*c3+c5*c4-c4*c3=0`); end if; if c4=0 then ERROR(`RK_45: c4=0`); end if; if c3=c4 then ERROR(`RK_45: c3=c4`); end if; if c3=c5 then ERROR(`RK_45: c3=c5`); end if; if c4=c5 then ERROR(`RK_45: c4=c5`); end if; a[2,1] := c2; a[3,2] := a32; a[3,1] := c3-a32; a[4,2] := a42; a[6,2] := b2; # Stage 1: b3, b4, b5 (è b1): a[6,3] := -1/12*(12*b2*c2^3+4*c5-3-12*c5*b2*c2^2-6*c5*c4-12*c4*b2*c2^2+12*c4*c5*b2*c2+4*c4)/c3/d; a[6,4] := 1/12*(6*c5*c3-4*c3-12*c3*c5*b2*c2+12*c3*b2*c2^2+3-12*b2*c2^3-4*c5+12*c5*b2*c2^2)/(c3-c4)/c4/(c5-c4); a[6,4] := factor(a[6,4]); if a[6,4]=0 then ERROR(`RK_45: a[6,4]=0`); end if; a[6,5] := -1/12/c5*(6*c4*c3+12*c3*b2*c2^2-12*c3*c4*b2*c2-4*c3-12*b2*c2^3+3+12*c4*b2*c2^2-4*c4)/(c5-c4)/(c3-c5); a[6,5] := factor(a[6,5]); if a[6,5]=0 then ERROR(`RK_45: a[6,5]=0`); end if; a[6,1] := 1-a[6,2]-a[6,3]-a[6,4]-a[6,5]; # Stage 2: a43 a[4,3] := -1/24*(-24*a[6,3]*a[3,2]*c2*c3 - 24*a[6,4]*c4*a[4,2]*c2 - 4*c5+24*c5*a[6,3]*a[3,2]*c2 +24*c5*a[6,4]*a[4,2]*c2 + 3) /c3/a[6,4]/(c5-c4); a[4,1] := c4-a[4,2]-a[4,3]; # rest only a[5,*] # Stage 3: a5* d1 := factor(-a[3,2]*c2^2*c4 + a[3,2]*c2*c4^2 + a[4,2]*c2^2*c3 + a[4,3]*c3^2*c2 -a[4,2]*c2*c3^2-a[4,3]); if d1=0 then ERROR(`RK_45: illegal denom. for a[5,4]=0`); end if; a[5,4] := -1/24*(4*a[3,2]*c2^2 + c3^2 - c2*c3 - 2*a[3,2]*c2)/ (-a[3,2]*c2^2*c4 + a[3,2]*c2*c4^2 + a[4,2]*c2^2*c3 + a[4,3]*c3^2*c2 -a[4,2]*c2*c3^2-a[4,3]*c3^3)/a[6,5]; a[5,3] := -1/24*(24*a[6,4]*a[4,3]*a[3,2]*c2 + 24*a[6,5]*a[5,4]*a[4,3]*c3 +24*a[6,5]*a[5,4]*a[4,2]*c2-1)/a[3,2]/c2/a[6,5]; a[5,2] := -1/6*(6*a[6,3]*a[3,2]*c2 + 6*a[6,4]*a[4,2]*c2 + 6*a[6,4]*a[4,3]*c3 +6*a[6,5]*a[5,3]*c3 + 6*a[6,5]*a[5,4]*c4-1)/a[6,5]/c2; a[5,1] := c5 - a[5,4] - a[5,3] - a[5,2]; simplify(a); end proc; # --------------------------------------------------------------------- RK_56 := proc(c2, c3, c4, c5, a43)::Matrix; local i,j, a, c6, b3,b4,b5,b6, c2_, c3_, c4_, c5_, a43_; c2_ := simplify(c2); c3_ := simplify(c3); c4_ := simplify(c4); c5_ := simplify(c5); a43_:= simplify(a43); if c2_= 0 then ERROR(`RK_56: c2= 0`); end if; if c3_= 0 then ERROR(`RK_56: c3= 0`); end if; if c3_= 1 then ERROR(`RK_56: c3= 1`); end if; if c4_= 0 then ERROR(`RK_56: c4= 0`); end if; if c4_= 1 then ERROR(`RK_56: c4= 1`); end if; if c5_= 0 then ERROR(`RK_56: c5= 0`); end if; if c5_= 1 then ERROR(`RK_56: c5= 1`); end if; # if c2_=c3_ then print(c2,c3);ERROR(`RK_56: c2=c3`); end if; # if c2_=c4_ then print(c2,c4);ERROR(`RK_56: c2=c4`); end if; # if c2_=c5_ then print(c2,c5);ERROR(`RK_56: c2=c5`); end if; if c3_=c4_ then print(c3,c4);ERROR(`RK_56: c3=c4`); end if; if c3_=c5_ then print(c3,c5);ERROR(`RK_56: c3=c5`); end if; if c4_=c5_ then print(c4,c5);ERROR(`RK_56: c4=c5`); end if; a := Matrix(7,7); a[7,2]:=0; a[4,3]:=a43_; a[2,1] := c2_; c6 := 1; b3:=(10*c4_*c5_-5*c4_-5*c5_+3)/c3_/(1-c3_)/(c3_-c5_)/(c3_-c4_)/60; b4:=(10*c5_*c3_-5*c3_-5*c5_+3)/(-1+c4_)/(c3_-c4_)/c4_/(c4_-c5_)/60; b5:=(3-5*c3_-5*c4_+10*c4_*c3_)/c5_/(1-c5_)/(c4_-c5_)/(c3_-c5_)/60; b6:=(30*c4_*c5_*c3_-20*c4_*c3_+15*c3_-20*c5_*c3_+15*c4_-20*c4_*c5_-12+15*c5_)/ (c3_-1)/(c5_-1)/(c4_-1)/60; b6:=factor(b6); if b6=0 then ERROR(`RK_56: b6=0`); end if; a[7,3]:=b3; a[7,4]:=b4; a[7,5]:=b5; a[7,6]:=b6; a[7,1]:=1-b3-b4-b5-b6; a[6,5] := factor(b5*(1-c5_)/b6); a[6,3] := (-12*c4_^3*b4+16*c4_^4*b4-4*b3*c3_^4-2*c4_*b3*c3_^4- 16*c4_*b5*c5_^4-4*c5_*c4_*b6+12*c4_^3*b4*c5_+ 12*a[4,3]*c3_^2*b4*c4_*c5_-4*c4_*b3*c3_^3*c5_-12*b4*c4_^2*c3_*c5_+ 4*c4_*b6-12*a[4,3]*c3_^2*b4*c4_^2+28*b5*c5_^3*c4_-16*b4*c4_^4*c5_+ 4*c5_*b3*c3_^4+4*c4_*b3*c3_^3+12*b4*c4_^2*c3_-12*b4*a[4,3]*c3_^3*c5_+ 12*a[4,3]*c3_^3*b4*c4_-2*c3_*b5*c5_^3*c4_+16*c5_*c3_*c4_^3*b4- 28*c3_*b5*c5_^3+17*c3_*b5*c5_^4-c4_^4*b4*c3_+4*c5_*c3_*b6- 2*c3_*c4_*b6-16*c4_^3*b4*c3_+12*b5*c5_^2*c3_-12*b5*c5_^2*c4_- 3*c3_*b6+b3*c3_^5)/(c5_-1)/(-c4_+c3_)/c3_^2/b6/12; a[5,4] := (-b5*c5_^4+2*c3_*b5*c5_^3-b6-c4_^4*b4+b3*c3_^4+ 2*c4_^3*b4*c3_+2*c3_*b6)/c4_/(c5_-1)/b5/(c4_-c3_)/12; a[5,3] := (b5*c5_^4-2*b5*c5_^3*c4_+b6+b3*c3_^4-12*b4*a[4,3]*c3_^2- 2*c4_*b3*c3_^3+12*a[4,3]*c3_^2*b4*c4_+12*b4*c4_*a[4,3]*c3_- c4_^4*b4-2*c4_*b6-12*b4*c4_^2*a[4,3]*c3_)/ (c5_-1)/(c4_-c3_)/b5/c3_/12; a[6,4] := factor((b4*(1-c4_)-b5*a[5,4])/b6); a[6,2]:= (c6^2/2 - a[6,3]*c3 - a[6,4]*c4 - a[6,5]*c5)/c2; a[5,2]:= (c5^2/2 - a[5,3]*c3 - a[5,4]*c4)/c2; a[4,2]:= (c4^2/2 - a[4,3]*c3)/c2; a[3,2]:= (c3^2/2)/c2; #for i from 3 to n do # A[i,2] := (c[i]^2/2 - add(A[i,j]*c[j], j=3..i-1) )/c[2]; #end do; a[3,1] := c3_ - a[3,2]; a[4,1] := c4_ - a[4,2] - a[4,3]; a[5,1] := c5_ - a[5,2] - a[5,3] - a[5,4]; a[6,1] := 1 - a[6,2] - a[6,3] - a[6,4] - a[6,5]; simplify(a); end proc; # --------------------------------------------------------------------- RK_56a := proc(c2, c3, c4, c5, b2, a43)::Matrix; local i,j, a, c6, b3,b4,b5,b6, den, u3,v3, u4,v4, u5,v5, A_,B_,C_, D_, d3,d4,d5,d6, t1,t2, sol, c2_, c3_, c4_, c5_, b2_, a43_; c2_ := simplify(c2); c3_ := simplify(c3); c4_ := simplify(c4); c5_ := simplify(c5); b2_ := simplify(b2); a43_:= simplify(a43); if c2_= 0 then ERROR(`RK_56a: c2= 0`); end if; if c2_=c3_ then ERROR(`RK_56a: c2=c3`); end if; if c2_=c4_ then ERROR(`RK_56a: c2=c4`); end if; if c2_=c5_ then ERROR(`RK_56a: c2=c5`); end if; if c3_=c4_ then ERROR(`RK_56a: c3=c4`); end if; if c3_=c5_ then ERROR(`RK_56a: c3=c5`); end if; if c4_=c5_ then ERROR(`RK_56a: c4=c5`); end if; if c5_= 1 then ERROR(`RK_56a: c5= 1`); end if; a := Matrix(7,7); a[7,2]:=b2_; a[4,3]:=a43_; a[2,1] := c2_; c6 := 1; sol := solve( { b2_*c2_ + b3*c3_ + b4*c4_ + b5*c5_ + b6*c6 - 1/2, b2_*c2_^2 + b3*c3_^2 + b4*c4_^2 + b5*c5_^2 + b6*c6^2 - 1/3, b2_*c2_^3 + b3*c3_^3 + b4*c4_^3 + b5*c5_^3 + b6*c6^3 - 1/4, b2_*c2_^4 + b3*c3_^4 + b4*c4_^4 + b5*c5_^4 + b6*c6^4 - 1/5}, { b3,b4,b5,b6} ); assign(sol); b3:=simplify(b3); b4:=simplify(b4); b5:=simplify(b5); b6:=simplify(b6); a[7,3]:=b3; a[7,4]:=b4; a[7,5]:=b5; a[7,6]:=b6; a[7,1]:=1-b2_-b3-b4-b5-b6; # find d3,d4,d5,d6 from equations: # b3*d3 + b4*d4 + b5*d5 + b6*d6 = 1/6; # b3*c3_*d3 + b4*c4_*d4 + b5*c5_*d5 + b6*d6 = 1/8; # b3*c3_^2*d3 + b4*c4_^2*d4 + b5*c5_^2*d5 + b6*d6 = 1/10; # b3*d3^2 + b4*d4^2 + b5*d5^2 + b6*d6^2 = 1/20; if b3=0 then ERROR(`RK_56a: b3=0`); end if; if b5=0 then ERROR(`RK_56a: b5=0`); end if; if b6=0 then ERROR(`RK_56a: b6=0`); end if; den := (c4_-c3_)*(c5_-c3_)*b3; u3 := (c4_*c5_/6 - c4_/8 - c5_/8 + 1/10)/den; v3 := -b6*(c5_-1)*(c4_-1)/den; den := (c4_-c3_)*(c4_-c5_)*b4; u4 := (c3_*c5_/6 - c5_/8 - c3_/8 + 1/10)/den; v4 := -b6*(c5_-1)*(c3_-1)/den; den := (c5_-c3_)*(c5_-c4_)*b5; u5 := (c3_*c4_/6 - c4_/8 - c3_/8 + 1/10)/den; v5 := -b6*(c3_-1)*(c4_-1)/den; # Now stays only one square equation: # b3*d3^2 + b4*d4^2 + b5*d5^2 + b6*d6^2 = 1/20; # or # d6^2*(b3*v3^2 + b4*v4^2 + b5*v5^2 + b6 ) + # d6*2*(b3*u3*v3 + b4*u4*v4 + b5*u5*v5) + # (b3*u3^2 + b4*u4^2 + b5*u5^2 - 1/20) # A*y^2 + B*y + C=0 A_ := simplify(b3*v3*v3 + b4*v4*v4 + b5*v5*v5 + b6); B_ := simplify(2*(b3*u3*v3 + b4*u4*v4 + b5*u5*v5)); C_ := simplify(b3*u3*u3 + b4*u4*u4 + b5*u5*u5 - 1/20); if A_=0 then ERROR(`RK_56a: A_=0`); end if; # Diskriminant: D_ := simplify(B_*B_ - 4*A_*C_); # if D_<0 then ERROR(`RK_56a: Discr_=0`); end if; #print(777, D_); #if D_<0 then a[1,1]:=999; return a; end if; D_ := sqrt(D_); d6 := simplify((-B_+D_)/(2*A_)); d3 := u3 + v3*d6; d4 := u4 + v4*d6; d5 := u5 + v5*d6; a[3,2] := simplify(d3/c2_); a[3,1] := simplify(c3_-a[3,2]); a[4,2] := simplify((d4-a43_*c3_)/c2_); a[4,1] := simplify(c4_ - a[4,2] - a43_); a[6,5] := simplify(b5*(1-c5_)/b6); if d4 =0 then ERROR(`RK_set56a: d4=0` ); end if; if a[6,5]=0 then ERROR(`RK_set56a: a[6,5]=0`); end if; t2 := simplify(b6*a[6,5]*(c3_*c3_*d4-c2_*c3_*d4-d3*c4_*c4_+c2_*d3*c4_)); if t2=0 then ERROR(`RK_56a: c3*c3*d4-c2*c3*d4-d3*c4*c4+c2*d3*c4=0`); end if; t1:=-2*d4+c4_*(c4_-c2_) +120*( c2_*c3_*c4_*d4*b4*a43_ -c2_*c4_*c4_*b4*d3*a43_ +c2_*d4*d5*b6*a[6,5] +c3_*c3_*b4*d4*a43_ -c3_*c3_*c4_*b4*d4*a43_ +c2_*c4_*b4*d3*a43_ -c2_*c3_*b3*d3*d4 +c2_*b3*d3*d4 +c2_*b4*d4*d4 -c2_*c4_*b4*d4*d4 -c2_*c3_*b4*d4*a43_ -c4_*c4_*b4*d3*a43_ +c4_*c4_*c4_*b4*d3*a43_ ); a[5,3] := simplify(-t1/(t2*120)); t2 := 120*b6*a[6,5]*d4; t1 := 1 -120*b4*a43_*d3 +120*b4*a43_*d3*c4_ -120*b6*a[6,5]*a[5,3]*d3; a[5,4] := simplify(t1/t2); a[5,2] := simplify((d5-a[5,3]*c3_-a[5,4]*c4_)/c2_); a[5,1] := c5_ - a[5,2] - a[5,3] - a[5,4]; a[6,4] := simplify(-(b5*a[5,4]-b4+b4*c4_)/b6); a[6,3] := simplify(-(b4*a[4,3]+b5*a[5,3]-b3+b3*c3_)/b6); a[6,2] := simplify((d6-a[6,3]*c3_-a[6,4]*c4_-a[6,5]*c5_)/c2_); a[6,1] := 1 - a[6,2] - a[6,3] - a[6,4] - a[6,5]; simplify(a); end proc; # # a6 := RK_56a(1/3, 1/2, 2/3, 3/4, 9/25, 1/3); # # [ 0 0 0 0 0 0 0] # [ ] # [1/3 0 0 0 0 0 0] # [ ] # [-5/8 9/8 0 0 0 0 0] # [ ] # [-2/3 1 1/3 0 0 0 0] # [ ] # [-179 75 25 15 ] # [---- --- -- --- 0 0 0] # [448 112 56 448 ] # [ ] # [ 112 ] # [7/6 7/10 5/3 -10 --- 0 0] # [ 15 ] # [ ] # [ -9 224 ] # [1/9 9/25 2/5 -- --- 1/30 0] # [ 10 225 ] # --------------------------------------------------------------------- # Discriminant in RK_56a: RK_56aD := proc(c2, c3, c4, c5, b2, a43) local i,j, a, c6, b3,b4,b5,b6, den, u3,v3, u4,v4, u5,v5, A_,B_,C_, D_, d3,d4,d5,d6, t1,t2, sol, c2_, c3_, c4_, c5_, b2_, a43_; c2_ := simplify(c2); c3_ := simplify(c3); c4_ := simplify(c4); c5_ := simplify(c5); b2_ := simplify(b2); a43_:= simplify(a43); if c2_= 0 then ERROR(`RK_56a: c2= 0`); end if; if c2_=c3_ then ERROR(`RK_56a: c2=c3`); end if; if c2_=c4_ then ERROR(`RK_56a: c2=c4`); end if; if c2_=c5_ then ERROR(`RK_56a: c2=c5`); end if; if c3_=c4_ then ERROR(`RK_56a: c3=c4`); end if; if c3_=c5_ then ERROR(`RK_56a: c3=c5`); end if; if c4_=c5_ then ERROR(`RK_56a: c4=c5`); end if; if c5_= 1 then ERROR(`RK_56a: c5= 1`); end if; a := Matrix(1,1); # a[7,2]:=b2_; a[4,3]:=a43_; a[2,1] := c2_; c6 := 1; sol := solve( { b2_*c2_ + b3*c3_ + b4*c4_ + b5*c5_ + b6*c6 - 1/2, b2_*c2_^2 + b3*c3_^2 + b4*c4_^2 + b5*c5_^2 + b6*c6^2 - 1/3, b2_*c2_^3 + b3*c3_^3 + b4*c4_^3 + b5*c5_^3 + b6*c6^3 - 1/4, b2_*c2_^4 + b3*c3_^4 + b4*c4_^4 + b5*c5_^4 + b6*c6^4 - 1/5}, { b3,b4,b5,b6} ); assign(sol); b3:=simplify(b3); b4:=simplify(b4); b5:=simplify(b5); b6:=simplify(b6); # a[7,3]:=b3; a[7,4]:=b4; a[7,5]:=b5; a[7,6]:=b6; # a[7,1]:=1-b2_-b3-b4-b5-b6; # find d3,d4,d5,d6 from equations: # b3*d3 + b4*d4 + b5*d5 + b6*d6 = 1/6; # b3*c3_*d3 + b4*c4_*d4 + b5*c5_*d5 + b6*d6 = 1/8; # b3*c3_^2*d3 + b4*c4_^2*d4 + b5*c5_^2*d5 + b6*d6 = 1/10; # b3*d3^2 + b4*d4^2 + b5*d5^2 + b6*d6^2 = 1/20; if b3=0 then ERROR(`RK_56a: b3=0`); end if; if b5=0 then ERROR(`RK_56a: b5=0`); end if; if b6=0 then ERROR(`RK_56a: b6=0`); end if; den := (c4_-c3_)*(c5_-c3_)*b3; u3 := (c4_*c5_/6 - c4_/8 - c5_/8 + 1/10)/den; v3 := -b6*(c5_-1)*(c4_-1)/den; den := (c4_-c3_)*(c4_-c5_)*b4; u4 := (c3_*c5_/6 - c5_/8 - c3_/8 + 1/10)/den; v4 := -b6*(c5_-1)*(c3_-1)/den; den := (c5_-c3_)*(c5_-c4_)*b5; u5 := (c3_*c4_/6 - c4_/8 - c3_/8 + 1/10)/den; v5 := -b6*(c3_-1)*(c4_-1)/den; # Now stays only one square equation: # b3*d3^2 + b4*d4^2 + b5*d5^2 + b6*d6^2 = 1/20; # or # d6^2*(b3*v3^2 + b4*v4^2 + b5*v5^2 + b6 ) + # d6*2*(b3*u3*v3 + b4*u4*v4 + b5*u5*v5) + # (b3*u3^2 + b4*u4^2 + b5*u5^2 - 1/20) # A*y^2 + B*y + C=0 A_ := simplify(b3*v3*v3 + b4*v4*v4 + b5*v5*v5 + b6); B_ := simplify(2*(b3*u3*v3 + b4*u4*v4 + b5*u5*v5)); C_ := simplify(b3*u3*u3 + b4*u4*u4 + b5*u5*u5 - 1/20); if A_=0 then ERROR(`RK_56a: A_=0`); end if; # Diskriminant: D_ := simplify(B_*B_ - 4*A_*C_); simplify(D_); end proc; # --------------------------------------------------------------------- # p=6 n=7 b2=0 Butcher solution RK_67 := proc(c2_, c3_, c5_, c6_)::Matrix; local n,n1, t1,t2, A, i,j, c4,c7,b2,b3,b4,b5,b6,b7, c2,c3,c5,c6, a32, a42,a43, a52,a53,a54, a62,a63,a64,a65, a72,a73,a74,a75,a76; n:=7:n1:=n+1: c2 := simplify(c2_); c3 := simplify(c3_); c5 := simplify(c5_); c6 := simplify(c6_); c7:=1: b2:=0: # ti - temporary variable t1 := 15*c3^2 - 10*c3 + 2: t2 := 5*c3^2 - 5*c3 + 1: c4 := c3/t1: a32 := 1/2*c3^2/c2: #print(790, a32); a42 := 1/2*(3*t1-2)*c3^2/t1^3/c2: a43 := -c3*(t1-1)/t1^3: a52 := 1/2*c3*c5*(5*c5^2-5*c5+1)/c2/t2: #print(778); a53 := -1/6*(150*c5*c3^3-145*c5*c3^2+40*c5*c3-4*c5+20*c3^2-c3-30*c3^3)* (c5-c3)*c5/c3^2/t2/(t1-1): a54 := -1/6*c5*(15*c5*c3^2-10*c5*c3+2*c5-c3)*(c5-c3)*t1^2/t2/(t1-1)/c3^2: a62 := 1/2*(-5*c6+5*c6^2+1)*c6*c3/c2/t2: # print(779); a65 := (3*c3-1)*(15*c6*c3^2-10*c6*c3+2*c6-c3)*(-c3+c6)*(-c6+c5)*c6/ (-9*c3+15*c5*c3+4-6*c5)/(c5-c3)/(15*c5*c3^2-10*c5*c3+2*c5-c3)/c5: b3 := -1/60*(-45*c5*c3^2+75*c5*c6*c3^2-45*c6*c3^2+30*c3^2+35*c5*c3-23*c3- 60*c5*c6*c3+35*c6*c3-6*c6+10*c5*c6-6*c5+4)/ c3^2/(15*c3^2-10*c3+1)/(c3-1)/(-c3+c6)/(c5-c3): b4 := 1/60*t1^5*(-5*c5*c3+3*c3-5*c6*c3+10*c5*c6*c3-5*c5*c6-2+3*c5+3*c6)/ c3^2/(15*c3^2-10*c3+1)/(5*c3-2)/(3*c3-1)/ (15*c6*c3^2-10*c6*c3+2*c6-c3)/(15*c5*c3^2-10*c5*c3+2*c5-c3): b5 := -1/60*(5*c3^2-5*c3+1)*(15*c6*c3-9*c3-6*c6+4)/c5/(c5-1)/ (15*c5*c3^2-10*c5*c3+2*c5-c3)/(c5-c3)/(-c6+c5): # print(780); b6 := 1/60*(5*c3^2-5*c3+1)*(-9*c3+15*c5*c3+4-6*c5)/(c6-1)/ (15*c6*c3^2-10*c6*c3+2*c6-c3)/(-c3+c6)/(-c6+c5)/c6: b7 := (1/2-b3*c3-b4*c4-b5*c5-b6*c6)/c7: a76 := factor((1-c6)*b6/b7): a75 := factor(((1-c5)*b5 - b6*a65)/b7): a72 := -(b6*a62+b5*a52+b4*a42+b3*a32)/b7: a64 := -(b6*a65*a54*a43*a32*c2+ b7*a75*a54*a43*a32*c2+ b7*a76*a65*a53*a32*c2+ b7*a76*a65*a54*a42*c2+ b7*a76*a65*a54*a43*c3 -1/720)/(b7*a76*a43*a32*c2): # print(782); # from 1st simplifying assumpion: a74 := factor(((1-c4)*b4 - b6*a64 - b5*a54)/b7): # from 2nd simplifying assumpion: a63 := factor(c6^2/2 - a62*c2 - a64*c4 - a65*c5)/c3: # from 1st simplifying assumpion: a73 := factor(((1-c3)*b3 - b6*a63 - b5*a53 - b4*a43)/b7): A := Matrix(n1,n1); A[2,1]:=c2: A[3,2]:=a32; A[4,2]:=a42; A[4,3]:=a43; A[5,2]:=a52; A[5,3]:=a53; A[5,4]:=a54; A[6,2]:=a62; A[6,3]:=a63; A[6,4]:=a64; A[6,5]:=a65; A[7,2]:=a72; A[7,3]:=a73; A[7,4]:=a74; A[7,5]:=a75; A[7,6]:=a76; A[8,2]:=b2; A[8,3]:=b3; A[8,4]:=b4; A[8,5]:=b5; A[8,6]:=b6; A[8,7]:=b7; A[3,1]:=c3-A[3,2]; A[4,1]:=c4-add( A[4,j], j=2..3); A[5,1]:=c5-add( A[5,j], j=2..4); A[6,1]:=c6-add( A[6,j], j=2..5); A[7,1]:=c7-add( A[7,j], j=2..6); A[8,1]:=1-add( A[8,j], j=2..7); simplify(A); end proc; # --------------------------------------------------------------------- # A := RK_67(1/2, 2/3, 5/6, 1/6); # # [ 0 0 0 0 0 0 0 0] # [ ] # [1/2 0 0 0 0 0 0 0] # [ ] # [2/9 4/9 0 0 0 0 0 0] # [ ] # [ -1 ] # [7/36 2/9 -- 0 0 0 0 0] # [ 12 ] # [ ] # [-35 -55 35 ] # [--- --- -- 15/8 0 0 0 0] # A := [144 36 48 ] # [ ] # [-1 -11 ] # [--- --- -1/8 1/2 1/10 0 0 0] # [360 36 ] # [ ] # [-41 22 43 -118 32 80 ] # [--- -- --- ---- --- -- 0 0] # [260 13 156 39 195 39 ] # [ ] # [13 11 11 13 ] # [--- 0 -- -- 4/25 4/25 --- 0] # [200 40 40 200 ] # --------------------------------------------------------------------- RK_67a_R :=(x, x1,x2,x3) -> (1 - 7/5*(x1+x2+x3) + 2*(x1*x2+x1*x3+x2*x3) - 3*x1*x2*x3) +x *(-14/5 + 4 *(x1+x2+x3) - 6*(x1*x2+x1*x3+x2*x3) + 10*x1*x2*x3) +x^2*(2 - 3 *(x1+x2+x3) + 5*(x1*x2+x1*x3+x2*x3) - 10*x1*x2*x3): RK_67a := proc(c3_, c5_, c6_)::Matrix; local a32, a42,a43, a52,a53,a54, a62,a63,a64,a65, a72,a73,a74,a75, a76, b1,b2,b3,b4,b5,b6,b7, c2,c3,c4,c5,c6,c7, d3,d4,d5,d6,d7, e3,e4,e5,e6,e7, eq09,eq14,eq29,eq32,eq33,eq34, s1,s2,s3,s4,b2a, A; c3 := factor(c3_); c5 := factor(c5_); c6 := factor(c6_); eq09 := 24*(b5*a54*a43*a32*c2+ b6*a64*a43*a32*c2+ b6*a65*a53*a32*c2+ b6*a65*a54*a42*c2+ b6*a65*a54*a43*c3+ b7*a74*a43*a32*c2+ b7*a75*a53*a32*c2+ b7*a75*a54*a42*c2+ b7*a75*a54*a43*c3+ b7*a76*a63*a32*c2+ b7*a76*a64*a42*c2+ b7*a76*a64*a43*c3+ b7*a76*a65*a52*c2+ b7*a76*a65*a53*c3+ b7*a76*a65*a54*c4) -b2*c2^4-b3*c3^4-b4*c4^4-b5*c5^4-b6*c6^4-b7; eq14 := 3*(b3*a32*c2^2*c3+ b4*c4*a42*c2^2+ b4*c4*a43*c3^2+ b5*c5*a52*c2^2+ b5*c5*a53*c3^2+ b5*c5*a54*c4^2+ b6*c6*a62*c2^2+ b6*c6*a63*c3^2+ b6*c6*a64*c4^2+ b6*c6*a65*c5^2+ b7*a72*c2^2+ b7*a73*c3^2+ b7*a74*c4^2+ b7*a75*c5^2+ b7*a76*c6^2) -b2*c2^4-b3*c3^4-b4*c4^4-b5*c5^4-b6*c6^4-b7; eq29 := 4*(b3*a32*c2^3*c3+ b4*c4*a42*c2^3+ b4*c4*a43*c3^3+ b5*c5*a52*c2^3+ b5*c5*a53*c3^3+ b5*c5*a54*c4^3+ b6*c6*a62*c2^3+ b6*c6*a63*c3^3+ b6*c6*a64*c4^3+ b6*c6*a65*c5^3+ b7*a72*c2^3+ b7*a73*c3^3+ b7*a74*c4^3+ b7*a75*c5^3+ b7*a76*c6^3) -b2*c2^4-b3*c3^4-b4*c4^4-b5*c5^4-b6*c6^4-b7+1/30; eq32 := 6*(b7*a72*c2^2*a74*c4+ b7*a72*c2^2*a75*c5+ b7*a72*c2^2*a76*c6+ b7*a73*c3^2*a72*c2+ b7*a73*c3^2*a74*c4+ b7*a72*c2^2*a73*c3+ b4*a42*c2^2*a43*c3+ b4*a43*c3^2*a42*c2+ b5*a52*c2^2*a53*c3+ b5*a52*c2^2*a54*c4+ b5*a53*c3^2*a52*c2+ b5*a53*c3^2*a54*c4+ b5*a54*c4^2*a52*c2+ b5*a54*c4^2*a53*c3+ b6*a62*c2^2*a63*c3+ b6*a62*c2^2*a64*c4+ b6*a62*c2^2*a65*c5+ b6*a63*c3^2*a62*c2+ b6*a63*c3^2*a64*c4+ b6*a63*c3^2*a65*c5+ b6*a64*c4^2*a62*c2+ b6*a64*c4^2*a63*c3+ b6*a64*c4^2*a65*c5+ b6*a65*c5^2*a62*c2+ b6*a65*c5^2*a63*c3+ b6*a65*c5^2*a64*c4+ b4*a42^2*c2^3+ b4*a43^2*c3^3+ b5*a52^2*c2^3+ b5*a53^2*c3^3+ b5*a54^2*c4^3+ b6*a62^2*c2^3+ b6*a63^2*c3^3+ b6*a64^2*c4^3+ b6*a65^2*c5^3+ b7*a72^2*c2^3+ b7*a73^2*c3^3+ b7*a74^2*c4^3+ b7*a75^2*c5^3+ b7*a76^2*c6^3+ b3*a32^2*c2^3+ b7*a73*c3^2*a75*c5+ b7*a73*c3^2*a76*c6+ b7*a74*c4^2*a72*c2+ b7*a74*c4^2*a73*c3+ b7*a74*c4^2*a75*c5+ b7*a74*c4^2*a76*c6+ b7*a75*c5^2*a72*c2+ b7*a75*c5^2*a73*c3+ b7*a75*c5^2*a74*c4+ b7*a75*c5^2*a76*c6+ b7*a76*c6^2*a72*c2+ b7*a76*c6^2*a73*c3+ b7*a76*c6^2*a74*c4+ b7*a76*c6^2*a75*c5) -b2*c2^4 -b3*c3^4 -b4*c4^4 -b5*c5^4 -b6*c6^4 -b7+1/30; eq33 := 6*(b7*a73*a32*c2+ b7*a74*a42*c2+ b7*a74*a43*c3+ b7*a75*a52*c2+ b7*a75*a53*c3+ b7*a75*a54*c4+ b7*a76*a62*c2+ b7*a76*a63*c3+ b7*a76*a64*c4+ b7*a76*a65*c5+ b5*c5^2*a53*a32*c2+ b5*c5^2*a54*a42*c2+ b5*c5^2*a54*a43*c3+ b6*c6^2*a63*a32*c2+ b6*c6^2*a64*a42*c2+ b6*c6^2*a64*a43*c3+ b6*c6^2*a65*a52*c2+ b6*c6^2*a65*a53*c3+ b6*c6^2*a65*a54*c4+ b4*a43*a32*c2*c4^2) -b2*c2^4 -b3*c3^4 -b4*c4^4 -b5*c5^4 -b6*c6^4 -b7+1/30; eq34 := 3*(b3*a32*c2^2*c3^2+ b4*c4^2*a42*c2^2+ b4*c4^2*a43*c3^2+ b5*c5^2*a52*c2^2+ b5*c5^2*a53*c3^2+ b5*c5^2*a54*c4^2+ b6*c6^2*a62*c2^2+ b6*c6^2*a63*c3^2+ b6*c6^2*a64*c4^2+ b6*c6^2*a65*c5^2+ b7*a72*c2^2+ b7*a73*c3^2+ b7*a74*c4^2+ b7*a75*c5^2+ b7*a76*c6^2) -b2*c2^4 -b3*c3^4 -b4*c4^4 -b5*c5^4 -b6*c6^4 -b7+1/30; # c4 :=(7 -99/5*c5+14*c5^2+c6*(-99/5+57*c5-41*c5^2)+c6^2*(14-41*c5+30*c5^2)) /(47/5-26*c5+18*c5^2+c6*(-26+72*c5-50*c5^2)+c6^2*(18-50*c5+35*c5^2)); # # # substitute c2(c3,c4,c5,c6): s1:=c3 + c4 + c5 + c6; s2:=c3*c4 + c3*c5 + c3*c6 + c4*c5 + c4*c6 + c5*c6; s3:=c3*c4*c5 + c3*c4*c6 + c3*c5*c6 + c4*c5*c6; s4:=c3*c4*c5*c6; c2:=factor((5-7*s1+10*s2-15*s3+25*s4)/(7-10*s1+15*s2-25*s3+50*s4)); b2a:=37500*RK_67a_R(c3,c4,c5,c6)* RK_67a_R(c6,c3,c4,c5)* RK_67a_R(c5,c6,c3,c4)* RK_67a_R(c4,c5,c6,c3); b2 := factor((7 - 10*s1 + 15*s2 - 25*s3 + 50*s4)^5/b2a); # b3 := 1/60*(5*c6*c5-10*c5*c4*c6-60*c6*b2*c2^4+60*c6*c5*b2*c2^3+ 60*c6*c4*b2*c2^3-60*c6*c5*b2*c2^2-60*c6*c4*b2*c2^2+ 60*c6*b2*c2^3+5*c4*c6-60*c6*c5*c4*b2*c2^2+ 60*c6*c5*b2*c2*c4-3*c6+2-3*c5-60*b2*c2^4+5*c4*c5+ 60*c5*b2*c2^3+60*c5*c4*b2*c2^3-60*c5*c4*b2*c2^2+60*b2*c2^5- 60*c4*b2*c2^4+60*c4*b2*c2^3-60*c5*b2*c2^4-3*c4)/ (c3-1)/c3/(-c3+c6)/(-c3+c5)/(-c3+c4); b3:=simplify(b3); # b4 := -1/60*(60*c6*c5*b2*c2^3-3*c6-60*c6*c3*b2*c2^2- 60*c6*c5*c3*b2*c2^2-10*c5*c3*c6+60*c6*b2*c2^3+ 60*c6*c3*b2*c2^3+60*c6*c5*c3*b2*c2-60*c6*c5*b2*c2^2- 60*c6*b2*c2^4+5*c3*c6+5*c6*c5+2-60*b2*c2^4-3*c3- 60*c3*b2*c2^4+60*b2*c2^5+60*c3*b2*c2^3-3*c5- 60*c5*c3*b2*c2^2+60*c5*b2*c2^3+60*c5*c3*b2*c2^3- 60*c5*b2*c2^4+5*c3*c5)/(c4-1)/(-c3+c4)/c4/(c5-c4)/(c6-c4); b4:=simplify(b4); # b5 := 1/60*(-3*c6-10*c3*c4*c6+5*c3*c6+5*c4*c6+60*c6*b2*c2^3- 60*c6*c4*c3*b2*c2^2+60*c6*c4*c3*b2*c2-60*c6*c3*b2*c2^2+ 60*c6*c4*b2*c2^3+60*c6*c3*b2*c2^3-60*c6*c4*b2*c2^2- 60*c6*b2*c2^4+2-60*c4*c3*b2*c2^2+5*c3*c4-60*b2*c2^4- 60*c3*b2*c2^4+60*b2*c2^5-60*c4*b2*c2^4+60*c4*b2*c2^3+ 60*c3*b2*c2^3-3*c4-3*c3+60*c3*c4*b2*c2^3)/ (c5-1)/(-c3+c5)/(c5-c4)/c5/(c6-c5); b5:=simplify(b5); # b6 := -1/60*(2-60*c4*c3*b2*c2^2-3*c5+5*c3*c4-60*b2*c2^4- 10*c3*c4*c5+5*c3*c5+5*c4*c5+60*c5*b2*c2^3-60*c3*b2*c2^4- 60*c5*c4*c3*b2*c2^2+60*c5*c4*c3*b2*c2 -60*c5*c3*b2*c2^2+60*c5*c4*b2*c2^3+60*c5*c3*b2*c2^3- 60*c5*c4*b2*c2^2+60*b2*c2^5 -60*c4*b2*c2^4+60*c4*b2*c2^3-60*c5*b2*c2^4+60*c3*b2*c2^3- 3*c4-3*c3+60*c3*c4*b2*c2^3) /c6/(c6-1)/(-c3+c6)/(c6-c4)/(c6-c5); b6:=simplify(b6); # e3 := 1/2*(c2-c5)*(c2-c4)*(-c2+c6)*b2*c2^2/b3/ (-c3+c6)/(c3-c5)/(c3-c4); e4 := -1/2*(c2-c5)*(c2-c3)*(-c2+c6)*b2*c2^2/(c3-c4)/ b4/(c6-c4)/(-c5+c4); e5 := 1/2*(c2-c4)*(c2-c3)*(-c2+c6)*b2*c2^2/(-c5+c4)/ (c3-c5)/b5/(c6-c5); e6 := 1/2*b2*c2^2*(c2-c5)*(c2-c4)*(c2-c3)/b6/ (c6-c5)/(c6-c4)/(-c3+c6); # e3 := simplify(e3); e4 := simplify(e4); e5 := simplify(e5); e6 := simplify(e6); # b7 := 1/12; b3:=simplify(b3); b4:=simplify(b4); b5:=simplify(b5); b6:=simplify(b6); # d3:=c3^2/2+e3; d4:=c4^2/2+e4; d5:=c5^2/2+e5; d6:=c6^2/2+e6; d7:=1/2; b1:=1/12; # a32:=factor(d3/c2); a42:=factor((d4-a43*c3)/c2); a52:=factor((d5-a53*c3-a54*c4)/c2); a62:=factor((d6-a63*c3-a64*c4-a65*c5)/c2); a72:=factor((d7-a73*c3-a74*c4-a75*c5-a76*c6)/c2); # a76 := factor((1-c6)*b6/b7); a75 := factor(((1-c5)*b5 - b6*a65)/b7); a74 := factor(((1-c4)*b4 - b6*a64 - b5*a54)/b7); a73 := factor(((1-c3)*b3 - b6*a63 - b5*a53 - b4*a43)/b7); # ================================================================================ eq09 := simplify(eval(eq09)); eq14 := simplify(eval(eq14)); eq29 := simplify(eval(eq29)); eq32 := simplify(eval(eq32)); eq33 := simplify(eval(eq33)); eq34 := simplify(eval(eq34)); a65 := solve(eq14,a65); eq09 := simplify(eval(eq09)); eq29 := simplify(eval(eq29)); eq32 := simplify(eval(eq32)); eq33 := simplify(eval(eq33)); eq34 := simplify(eval(eq34)); a53 := solve(eq34,a53); eq32 := simplify(eval(eq32)); a43 := solve(eq32,a43); eq09 := simplify(eval(eq09)); eq29 := simplify(eval(eq29)); eq33 := simplify(eval(eq33)); a64 := solve(eq29,a64); eq09 := simplify(eval(eq09)); eq33 := simplify(eval(eq33)); a54 := solve(eq33,a54); eq09 := simplify(eval(eq09)); a63 := solve(eq09,a63); A := Matrix(8,8); A[2,1]:=eval(c2 ): A[3,2]:=eval(a32); A[4,2]:=eval(a42); A[4,3]:=eval(a43); A[5,2]:=eval(a52); A[5,3]:=eval(a53); A[5,4]:=eval(a54); A[6,2]:=eval(a62); A[6,3]:=eval(a63); A[6,4]:=eval(a64); A[6,5]:=eval(a65); A[7,2]:=eval(a72); A[7,3]:=eval(a73); A[7,4]:=eval(a74); A[7,5]:=eval(a75); A[7,6]:=eval(a76); A[8,2]:=eval(b2 ); A[8,3]:=eval(b3 ); A[8,4]:=eval(b4 ); A[8,5]:=eval(b5 ); A[8,6]:=eval(b6); A[8,7]:=eval(b7); A[3,1]:=c3-A[3,2]; A[4,1]:=c4-add( A[4,j], j=2..3); A[5,1]:=c5-add( A[5,j], j=2..4); A[6,1]:=c6-add( A[6,j], j=2..5); A[7,1]:=1 -add( A[7,j], j=2..6); A[8,1]:=1-add( A[8,j], j=2..7); simplify(A); end: #---------------------------------------------------------------- # RK_67b(c2, a32, b3) ñ2=ñ3=ñ6 RK_67b := proc(c2_, a32_, b3_)::Matrix; local a32, a42,a43, a52,a53,a54, a62,a63,a64,a65, a72,a73,a74,a75, a76, b1,b2,b3,b4,b5,b6,b7, c2,c3,c4,c5,c6,c7, t1, t2, t3, t4, A, a, b, c; c2 := factor(c2_ ); a32:= factor(a32_); b3 := factor(b3_ ); c3 := c2; c6:=c2; t1 := (15*c2^2 - 10*c2 + 2); t2 := ( 5*c2^2 - 5*c2 + 1); c4 := c2/t1; c5 := (9*c2-4)/(5*c2-2)/3; b4 := 1/60*t1^5/(t1-1)/(3*c2-1)^2/ (45*c2^2-40*c2+8)/c2^2/(5*c2-2); b5 := 81/40*(5*c2-2)^5*(5*c2^2-5*c2+1)/(15*c2^2-15*c2+4)/ (3*c2-1)^2/(45*c2^2-40*c2+8)/(-4+9*c2); b7 := 1/120*(675*c2^4-1305*c2^3+865*c2^2-240*c2+24)/ (c2-1)/(3*c2-1)^2/(5*c2-2); a42 := 1/2*c4^2*(a32+c2-c4)/a32/c2; a43 := -1/2*c4^2*(c2-c4)/a32/c2; a54 := 1/162*(9*c2-4)*(45*c2^2-40*c2+8)*t1^2* (3*c2-1)*(3*t2+1)/c2^2/(t1-1)/t2/(5*c2-2)^3; a52 := -1/2*(c5^3-c5^2*a32-c2*c5^2-3*a54*c4^2+2*c2*a54*c4+2*a54*c4*a32)/a32/c2; a53 := 1/2*(-3*a54*c4^2+2*c2*a54*c4-c2*c5^2+c5^3)/a32/c2; # t3, t4 - temp. var.: t3 := 120*b3*c2*(c2-1)*(t1-1)*(3*t2+1)*a32+3*t1-2; t4 := 240*b3*(c2-1)*(3*t2+1)*(t1-1)*a32^2 +3*t1-2; b6 := -1/60*t3^2/c2^2/(c2-1)/(3*t2+1)/(t1-1)/t4; b2 := -b3*(3*t1-2)*(2*a32-c2)^2/c2^2/t4; a62 := (t4/2*c2^2 + 1/2*(15*c2^2-20*c2+6)*(30*c2^2-15*c2+2)*(-2*a32+c2))*t4/2/a32/t3^2; a63 := 1/40*(5*c2-2)*(3*c2-1)*t2/c2/a32/(c2-1)/(3*t2+1)/(t1-1)/b6; a64 := 1/20*t2*(5*c2-2)*t1^2/(45*c2^2-40*c2+8)/c2^2/b6/(c2-1)/(t1-1); a65 := -9/20*t2*(5*c2-2)^3/b6/(c2-1)/(9*c2-4)/(45*c2^2-40*c2+8)/(3*t2+1); a76 := -b6*(-1+c2)/b7; # from eq43 a75 := -(b6*a65-b5+b5*c5)/b7; a74 := -(b4*c4+a64*b6+b5*a54-b4)/b7; a72 := -1/2*(-2*a76*c2*a65*c5-a32*c2+2*c2*a75*c5*a32+2*c2*a74*c4*a32- 2*a76*a62*c2*a32+2*a76*c2*a32^2-2*a76*c2*a64*c4-a32^2+a32- 3*a75*c5^2*a32+3*a76*a65*c5^2+2*a74*c4*a32^2-3*a74*c4^2*a32- 2*a76*a65*c5*a32+2*a75*c5*a32^2+3*a76*a64*c4^2-2*a76*a64*c4*a32) /a32^2/c2; a73 := 1/2*(-2*a76*c2*a65*c5-3*a74*c4^2*a32+3*a76*a65*c5^2-3*a75*c5^2*a32- 2*a76*c2*a64*c4+2*c2*a75*c5*a32+2*c2*a74*c4*a32-2*a76*a62*c2*a32- a32*c2+a32-2*a76*a65*c5*a32-2*a76*a64*c4*a32+3*a76*a64*c4^2) /a32^2/c2; A := Matrix(8,8); A[2,1]:=eval(c2 ): A[3,2]:=eval(a32); A[4,2]:=eval(a42); A[4,3]:=eval(a43); A[5,2]:=eval(a52); A[5,3]:=eval(a53); A[5,4]:=eval(a54); A[6,2]:=eval(a62); A[6,3]:=eval(a63); A[6,4]:=eval(a64); A[6,5]:=eval(a65); A[7,2]:=eval(a72); A[7,3]:=eval(a73); A[7,4]:=eval(a74); A[7,5]:=eval(a75); A[7,6]:=eval(a76); A[8,2]:=eval(b2 ); A[8,3]:=eval(b3 ); A[8,4]:=eval(b4 ); A[8,5]:=eval(b5 ); A[8,6]:=eval(b6); A[8,7]:=eval(b7); A[3,1]:=c3-A[3,2]; A[4,1]:=c4-add( A[4,j], j=2..3); A[5,1]:=c5-add( A[5,j], j=2..4); A[6,1]:=c6-add( A[6,j], j=2..5); A[7,1]:=1 -add( A[7,j], j=2..6); A[8,1]:=1-add( A[8,j], j=2..7); simplify(A); end: #---------------------------------------------------------------- RK_79 := proc(c4,c5)::Matrix; local s,c,b,a,A,uu,vv, t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12, eq13,eq28,eq65,eq75,eq81, sol,i,j; s:=9; c[4]:=c4: c[5]:=c5: c[1]:=0: c[s]:=1: c[s+1]:=1: c[2] :=c[4]/3: c[3] :=2/3*c[4]: a[3,2]:=2/3*c[4]: a[4,3]:=3/4*c[4]: uu := c[4]+c[5]: vv := c[4]*c[5]: t1 := 3-12*uu+24*vv+14*uu^2-70*uu*vv+105*vv^2; t2 := 5*c[4]*c[5]-2*c[4]-2*c[5]+1; t3 := 21*c[4]*c[5]-7*c[4]-7*c[5]+3; t4 := c[4]+c[5]-12*c[4]*c[5]+7*c[4]^2*c[5]+7*c[4]*c[5]^2; t5 := 105*c[4]^3*c[5]^2 -70*c[4]^3*c[5] -70*c[4]^2*c[5]^2 +45*c[4]^2*c[5] +14*c[4]^3 -12*c[4]^2 +7*c[4]*c[5]^2 +2*c[4]-c[5]; t6 := 105*c[4]^2*c[5]^3 -70*c[5]^3*c[4] -70*c[4]^2*c[5]^2 +45*c[4]*c[5]^2 +14*c[5]^3 -12*c[5]^2 +7*c[4]^2*c[5] +2*c[5]-c[4]; t7 := 315*c[4]^2*c[5]^2-280*c[4]^2*c[5]+56*c[4]^2-60*c[4]+15+274*c[4]*c[5]-60*c[5]-280*c[4]*c[5]^2+56*c[5]^2; t8 := 105*c[4]^2*c[5] -42*c[4]^2 -105*c[4]*c[5] +49*c[4] +28*c[5] -15; t9 := 105*c[4]*c[5]^2 -105*c[4]*c[5] +28*c[4] -42*c[5]^2 +49*c[5] -15; t10:= 63*c[4]*c[5]-28*c[4]-28*c[5]+15; t11:= 2-8*c[4]+7*c[4]^2-8*c[5]+37*c[4]*c[5]-35*c[4]^2*c[5]+7*c[5]^2-35*c[4]*c[5]^2+35*c[4]^2*c[5]^2; t12:= 9-84*c[4]-84*c[5]+287*c[4]^2+287*c[5]^2+862*c[4]*c[5]+13522*c[4]^2*c[5]^2-3246*c[4]^2*c[5] -3246*c[4]*c[5]^2-22050*c[4]^3*c[5]^2+4914*c[4]^3*c[5]-22050*c[4]^2*c[5]^3+4914*c[5]^3*c[4] -406*c[4]^3-406*c[5]^3-20580*c[5]^4*c[4]^3+11515*c[5]^4*c[4]^2-2450*c[5]^4*c[4] -20580*c[4]^4*c[5]^3+11025*c[4]^4*c[5]^4+38150*c[4]^3*c[5]^3+11515*c[4]^4*c[5]^2 -2450*c[4]^4*c[5]+196*c[4]^4+196*c[5]^4; c[8] :=0: b[2] :=0: b[3] :=0: a[4,2]:=0: a[5,2]:=0: a[6,2]:=0: a[7,2]:=0: a[8,2]:=0: a[9,2]:=0: a[8,3]:=0: b[9] := 1-b[4]-b[5]-b[6]-b[7]-b[8]: # -------------------------------------------- b[6] := t1^6/t2/t3^2/t4/t5/t6/t7/60; b[7] := 4084101*t11*t2^6/t9/t10/t3^2/t8/t7/40; b[8] := -t12/t4/c[5]/t10/c[4]/60; c[6] := (uu-12*vv+7*uu*vv)/t1; c[6] := factor(c[6]); c[7] := t10/t2/21; b[4] := (3*c[5]+6*c[5]*b[7]*c[7]-6*b[8]*c[5]+6*c[5]*b[6]*c[6]-6*b[7]*c[5]-6*b[6]*c[5]-6*b[6]*c[6]^2 -6*b[7]*c[7]^2-1+6*b[7]*c[7]+6*b[6]*c[6])/(c[4]-1)/(-c[5]+c[4])/6; b[5] := (-1+6*b[7]*c[7]-6*b[6]*c[6]^2-6*b[7]*c[7]^2+3*c[4]-6*c[4]*b[8]+6*c[4]*b[6]*c[6] +6*c[4]*b[7]*c[7]-6*c[4]*b[7]+6*b[6]*c[6]-6*c[4]*b[6])/(c[5]^2-c[4]*c[5]+c[4]-c[5])/6; a[5,4]:= c[5]^2*(c[5]-c[4]) /c[4]^2; a[5,3]:= 3*c[5]^2*(3*c[4]-2*c[5])/4/c[4]^2; a[6,5]:= (7*c[4]^2-6*c[4]+1)*t5*t6*t4/c[5]/(c[5]-c[4])/t1^4; a[7,6]:= t10*t9*t8*t1^3*t3*t7/t4/t2^4/t5/t6/t11/1166886; a[8,6]:= 3*c[5]*t2*t1^3*t11*t10*c[4]/t7/t6/t5/t12; a[8,5]:= 3*c[4]*t3*t2*t4*t10*t11/t12/(c[4]-c[5])/t9/t6; a[8,7]:= -27783*c[4]*c[5]*t11*t2^4*t4/t12/t7/t9/t8; a[8,4]:= 3*t11*t2*t4*t3*t10*c[5]/t12/t5/(c[5]-c[4])/t8; # 1st simpl.assumption: a[9,8] := b[8]/b[9]; a[9,7] := -(b[8]*a[8,7]-b[7]+b[7]*c[7])/b[9]; a[9,6] := -(b[7]*a[7,6]+b[8]*a[8,6]-b[6]+b[6]*c[6])/b[9]; a[9,5] := -(b[6]*a[6,5]+b[7]*a[7,5]+b[8]*a[8,5]-b[5]+b[5]*c[5])/b[9]; a[9,4] := -(b[5]*c[5]^3-b[5]*c[5]^2*c[4]+b[6]*c[4]^2*a[6,4]+b[7]*c[4]^2*a[7,4] +b[8]*c[4]^2*a[8,4]-b[4]*c[4]^2+b[4]*c[4]^3)/c[4]^2/b[9]; a[9,3] := -1/4*(3*b[4]*c[4]^3+9*b[5]*c[5]^2*c[4]-6*b[5]*c[5]^3 +4*b[6]*c[4]^2*a[6,3]+4*b[7]*c[4]^2*a[7,3])/c[4]^2/b[9]; eq13:=-15*b[4]*c[4]*a[9,5]*c[5]^2-15*b[4]*c[4]*a[9,6]*c[6]^2-15*b[4]*c[4]*a[9,7]*c[7]^2 -20/3*b[5]*c[5]*a[9,3]*c[4]^2-15*b[5]*c[5]*a[9,4]*c[4]^2-15*b[5]*c[5]*a[9,6]*c[6]^2 -15*b[5]*c[5]*a[9,7]*c[7]^2-20/3*b[6]*c[6]*a[9,3]*c[4]^2-15*b[6]*c[6]*a[9,4]*c[4]^2 -15*b[6]*c[6]*a[9,5]*c[5]^2-15*b[6]*c[6]*a[9,7]*c[7]^2-20/3*b[7]*c[7]*a[9,3]*c[4]^2 -15*b[7]*c[7]*a[9,4]*c[4]^2-15*b[7]*c[7]*a[9,5]*c[5]^2-15*b[7]*c[7]*a[9,6]*c[6]^2 -15*b[4]*c[4]^3*a[9,4]-15*b[6]*c[6]^3*a[9,6]-15*b[7]*c[7]^3*a[9,7] +5*b[4]*c[4]+5*b[5]*c[5]+5*b[6]*c[6]+5*b[7]*c[7]-15*b[5]*c[5]^3*a[9,5] -20/3*b[4]*c[4]^3*a[9,3]+10/3*a[9,3]*c[4]^2+15/2*a[9,4]*c[4]^2+15/2*a[9,5]*c[5]^2 +15/2*a[9,6]*c[6]^2+15/2*a[9,7]*c[7]^2+20/3*b[5]*c[5]*c[4]^2*a[5,3] +15*b[5]*c[5]*c[4]^2*a[5,4]+20/3*b[6]*c[6]*a[6,3]*c[4]^2+15*b[6]*c[6]*a[6,4]*c[4]^2 +15*b[6]*c[6]*a[6,5]*c[5]^2+20/3*b[7]*c[7]*a[7,3]*c[4]^2+15*b[7]*c[7]*a[7,4]*c[4]^2 +15*b[7]*c[7]*a[7,5]*c[5]^2+15*b[7]*c[7]*a[7,6]*c[6]^2-5*b[5]*c[5]^4-5*b[6]*c[6]^4 -5*b[7]*c[7]^4-5/2; eq28:=6*b[4]*c[4]+6*b[5]*c[5]+6*b[6]*c[6]+6*b[7]*c[7]-24*b[4]*c[4]^4*a[9,4]-64/9*b[4]*c[4]^4*a[9,3] +12*a[9,4]*c[4]^3+32/9*a[9,3]*c[4]^3+12*a[9,5]*c[5]^3+12*a[9,6]*c[6]^3+12*a[9,7]*c[7]^3 -24*b[5]*c[5]*a[9,4]*c[4]^3-24*b[6]*c[6]*a[9,4]*c[4]^3-24*b[7]*c[7]*a[9,4]*c[4]^3 -24*b[4]*c[4]*a[9,5]*c[5]^3-24*b[4]*c[4]*a[9,6]*c[6]^3-24*b[4]*c[4]*a[9,7]*c[7]^3 -64/9*b[5]*c[5]*a[9,3]*c[4]^3-24*b[5]*c[5]*a[9,6]*c[6]^3-24*b[5]*c[5]*a[9,7]*c[7]^3 -64/9*b[6]*c[6]*a[9,3]*c[4]^3-24*b[6]*c[6]*a[9,5]*c[5]^3-24*b[6]*c[6]*a[9,7]*c[7]^3 -64/9*b[7]*c[7]*a[9,3]*c[4]^3-24*b[7]*c[7]*a[9,5]*c[5]^3-24*b[7]*c[7]*a[9,6]*c[6]^3 -24*b[5]*c[5]^4*a[9,5]-24*b[6]*c[6]^4*a[9,6]-24*b[7]*c[7]^4*a[9,7]-6*b[4]*c[4]^4 -6*b[5]*c[5]^4-6*b[6]*c[6]^4-6*b[7]*c[7]^4+16/3*b[4]*c[4]^5+24*b[5]*a[5,4]*c[4]^3*c[5] +24*b[6]*c[6]*c[4]^3*a[6,4]+24*b[7]*c[7]*a[7,4]*c[4]^3+64/9*b[5]*c[5]*c[4]^3*a[5,3] +64/9*b[6]*c[6]*a[6,3]*c[4]^3+24*b[6]*c[6]*a[6,5]*c[5]^3+64/9*b[7]*c[7]*a[7,3]*c[4]^3 +24*b[7]*c[7]*a[7,5]*c[5]^3+24*b[7]*c[7]*a[7,6]*c[6]^3-14/5; eq65:=7*b[4]*c[4]+7*b[5]*c[5]+7*b[6]*c[6]+7*b[7]*c[7]-35*b[4]*c[4]^5*a[9,4]+35/2*a[9,4]*c[4]^4 +280/81*a[9,3]*c[4]^4+35/2*a[9,5]*c[5]^4+35/2*a[9,6]*c[6]^4+35/2*a[9,7]*c[7]^4 -560/81*b[4]*c[4]^5*a[9,3]-35*b[4]*c[4]*a[9,5]*c[5]^4-35*b[4]*c[4]*a[9,6]*c[6]^4 -35*b[4]*c[4]*a[9,7]*c[7]^4-35*b[5]*c[5]*a[9,6]*c[6]^4-35*b[5]*c[5]*a[9,7]*c[7]^4 -35*b[6]*c[6]*a[9,5]*c[5]^4-35*b[6]*c[6]*a[9,7]*c[7]^4-35*b[7]*c[7]*a[9,5]*c[5]^4 -35*b[7]*c[7]*a[9,6]*c[6]^4-35*b[5]*c[5]*a[9,4]*c[4]^4-35*b[6]*c[6]*a[9,4]*c[4]^4 -35*b[7]*c[7]*a[9,4]*c[4]^4-560/81*b[5]*c[5]*a[9,3]*c[4]^4-560/81*b[6]*c[6]*a[9,3]*c[4]^4 -560/81*b[7]*c[7]*a[9,3]*c[4]^4-35*b[5]*c[5]^5*a[9,5]-35*b[6]*c[6]^5*a[9,6]-35*b[7]*c[7]^5*a[9,7] -49/27*b[4]*c[4]^6+35*b[5]*a[5,4]*c[4]^4*c[5]+35*b[6]*c[6]*c[4]^4*a[6,4]+35*b[7]*c[7]*a[7,4]*c[4]^4 +560/81*b[5]*c[5]*c[4]^4*a[5,3]+560/81*b[6]*c[6]*c[4]^4*a[6,3]+560/81*b[7]*c[7]*a[7,3]*c[4]^4 +35*b[6]*c[6]*a[6,5]*c[5]^4+35*b[7]*c[7]*a[7,5]*c[5]^4+35*b[7]*c[7]*a[7,6]*c[6]^4-7*b[5]*c[5]^6 -7*b[6]*c[6]^6-7*b[7]*c[7]^6-7/2; eq75:=7*b[4]*c[4]+7*b[5]*c[5]+7*b[6]*c[6]+7*b[7]*c[7]-28*b[4]*c[4]^4*a[9,4]-224/27*b[4]*c[4]^4*a[9,3] +14*a[9,4]*c[4]^3+112/27*a[9,3]*c[4]^3+14*a[9,5]*c[5]^3+14*a[9,6]*c[6]^3+14*a[9,7]*c[7]^3 -28*b[5]*c[5]*a[9,4]*c[4]^3-28*b[6]*c[6]*a[9,4]*c[4]^3-28*b[7]*c[7]*a[9,4]*c[4]^3 -28*b[4]*c[4]*a[9,5]*c[5]^3-28*b[4]*c[4]*a[9,6]*c[6]^3-28*b[4]*c[4]*a[9,7]*c[7]^3 -224/27*b[5]*c[5]*a[9,3]*c[4]^3-28*b[5]*c[5]*a[9,6]*c[6]^3-28*b[5]*c[5]*a[9,7]*c[7]^3 -224/27*b[6]*c[6]*a[9,3]*c[4]^3-28*b[6]*c[6]*a[9,5]*c[5]^3-28*b[6]*c[6]*a[9,7]*c[7]^3 -224/27*b[7]*c[7]*a[9,3]*c[4]^3-28*b[7]*c[7]*a[9,5]*c[5]^3-28*b[7]*c[7]*a[9,6]*c[6]^3 -28*b[5]*c[5]^4*a[9,5]-28*b[6]*c[6]^4*a[9,6]-28*b[7]*c[7]^4*a[9,7]-7/9*b[4]*c[4]^6 +28*b[5]*a[5,4]*c[4]^3*c[5]^2+28*b[6]*c[4]^3*c[6]^2*a[6,4]+28*b[7]*c[7]^2*a[7,4]*c[4]^3 +224/27*b[5]*c[4]^3*c[5]^2*a[5,3]+224/27*b[6]*c[6]^2*a[6,3]*c[4]^3+28*b[6]*c[6]^2*a[6,5]*c[5]^3 +224/27*b[7]*c[7]^2*a[7,3]*c[4]^3+28*b[7]*c[7]^2*a[7,5]*c[5]^3+28*b[7]*c[7]^2*a[7,6]*c[6]^3 -7*b[5]*c[5]^6-7*b[6]*c[6]^6-7*b[7]*c[7]^6-7/2; eq81:=-21*b[4]*c[4]*a[9,5]*c[5]^2-21*b[4]*c[4]*a[9,6]*c[6]^2-21*b[4]*c[4]*a[9,7]*c[7]^2 -28/3*b[5]*c[5]*a[9,3]*c[4]^2-21*b[5]*c[5]*a[9,4]*c[4]^2-21*b[5]*c[5]*a[9,6]*c[6]^2 -21*b[5]*c[5]*a[9,7]*c[7]^2-28/3*b[6]*c[6]*a[9,3]*c[4]^2-21*b[6]*c[6]*a[9,4]*c[4]^2 -21*b[6]*c[6]*a[9,5]*c[5]^2-21*b[6]*c[6]*a[9,7]*c[7]^2-28/3*b[7]*c[7]*a[9,3]*c[4]^2 -21*b[7]*c[7]*a[9,4]*c[4]^2-21*b[7]*c[7]*a[9,5]*c[5]^2-21*b[7]*c[7]*a[9,6]*c[6]^2 -21*b[4]*c[4]^3*a[9,4]-21*b[6]*c[6]^3*a[9,6]-21*b[7]*c[7]^3*a[9,7]+7*b[4]*c[4]+7*b[5]*c[5] +7*b[6]*c[6]+7*b[7]*c[7]-21*b[5]*c[5]^3*a[9,5]-28/3*b[4]*c[4]^3*a[9,3]+14/3*a[9,3]*c[4]^2 +21/2*a[9,4]*c[4]^2+21/2*a[9,5]*c[5]^2+21/2*a[9,6]*c[6]^2+21/2*a[9,7]*c[7]^2 +28/3*b[5]*c[5]^3*c[4]^2*a[5,3]+21*b[5]*c[5]^3*c[4]^2*a[5,4]+28/3*b[6]*c[6]^3*a[6,3]*c[4]^2 +21*b[6]*c[6]^3*a[6,4]*c[4]^2+21*b[6]*c[6]^3*a[6,5]*c[5]^2+28/3*b[7]*c[7]^3*a[7,3]*c[4]^2 +21*b[7]*c[7]^3*a[7,4]*c[4]^2+21*b[7]*c[7]^3*a[7,5]*c[5]^2+21*b[7]*c[7]^3*a[7,6]*c[6]^2 -7*b[5]*c[5]^6-7*b[6]*c[6]^6-7*b[7]*c[7]^6-7/2; eq13:=eval(eq13); eq28:=eval(eq28); eq65:=eval(eq65); eq75:=eval(eq75); eq81:=eval(eq81); sol:=solve({eq13,eq28,eq65,eq75,eq81},{a[6,4],a[6,3],a[7,5],a[7,4],a[7,3]}); assign(sol); A := Matrix(10,10); for i from 2 to s do for j from 1 to i-1 do A[i,j]:=eval(a[i,j]); end do; end do; for j from 1 to s do A[s+1,j]:=eval(b[j]); end do; for i from 2 to s+1 do A[i,1]:=eval(c[i] - add(A[i,j], j=2..i-1)); end do; simplify(A); end proc: #---------------------------------------------------------------- RK_8_11 := proc(a105,b8)::Matrix; local A; A := Matrix(12,12): A[ 2, 1] := 1/2; A[ 3, 1] := 1/4; A[ 3, 2] := 1/4; A[ 4, 1] := 1/7; A[ 4, 2] := -1/14-3/98*21^(1/2); A[ 4, 3] := 5/49*21^(1/2)+3/7; A[ 5, 1] := 11/84+1/84*21^(1/2); A[ 5, 2] := 0; A[ 5, 3] := 2/7+4/63*21^(1/2); A[ 5, 4] := -1/252*21^(1/2)+1/12; A[ 6, 1] := 5/48+1/48*21^(1/2); A[ 6, 2] := 0; A[ 6, 3] := 1/36*21^(1/2)+1/4; A[ 6, 4] := -77/120+7/180*21^(1/2); A[ 6, 5] := 63/80-7/80*21^(1/2); A[ 7, 1] := 5/21-1/42*21^(1/2); A[ 7, 2] := 0; A[ 7, 3] := -24/35*a[10, 5]-136/105-(12/245*a[10,5])*sqrt(21)+656/2205*sqrt(21); A[ 7, 4] := 7-(3/10*a[10, 5])*sqrt(21)-71/45*sqrt(21)+3/10*a[10, 5]; A[ 7, 5] := -3/10*a[10, 5]+(3/10*a[10,5])*sqrt(21)-43/6+169/105*sqrt(21); A[ 7, 6] := -277/735*sqrt(21)+181/105+(12/245*a[10,5])*sqrt(21)+24/35*a[10, 5]; A[ 8, 1] := -(1/7560)*(180*b[8]*sqrt(21)-49*sqrt(21)-1800*b[8]+343)/b[8]; A[ 8, 2] := 0; A[ 8, 3] := -(1/99225)*(4860*a[10, 5]*sqrt(21)*b[8]+3087*a[10,5]*sqrt(21)-29520*b[8]*sqrt(21) +68040*a[10,5]*b[8]+7546*sqrt(21)-33957*a[10, 5]+128520*b[8]-33271)/b[8]; A[ 8, 4] :=-(1/16200)*(4860*a[10, 5]*sqrt(21)*b[8]-3528*a[10,5]*sqrt(21)+25560*b[8]*sqrt(21)-4860*a[10,5]*b[8]-6713*sqrt(21)+12348*a[10, 5]-113400*b[8]+29645)/b[8]; A[ 8, 5] := (1/56700)*(17010*a[10, 5]*sqrt(21)*b[8]-12348*a[10,5]*sqrt(21)+91260*b[8]*sqrt(21)-17010*a[10,5]*b[8]-24353*sqrt(21)+43218*a[10, 5]-406350*b[8]+108045)/b[8]; A[ 8, 6] := (1/198450)*(9720*a[10, 5]*sqrt(21)*b[8]+6174*a[10,5]*sqrt(21)-74790*b[8]*sqrt(21)+136080*a[10,5]*b[8]+16807*sqrt(21)-67914*a[10, 5]+342090*b[8]-78547)/b[8]; A[ 8, 7] := 49/1620/b[8]; A[ 9, 1] := 1/32; A[ 9, 2] := 0; A[ 9, 3] := (1/8)*a[10, 5]*sqrt(21)-(1/8)*a[10,5]-(1/72)*sqrt(21)+1/72; A[ 9, 4] := -49/288-(7/32*a[10,5])*sqrt(21)+7/288*sqrt(21)+49/32*a[10, 5]; A[ 9, 5] := (7/32*a[10, 5])*sqrt(21)-35/576*sqrt(21)-49/32*a[10,5]+21/64; A[ 9, 6] := -(1/8)*a[10, 5]*sqrt(21)+(1/8)*a[10,5]+(1/72)*sqrt(21)+5/36; A[ 9, 7] := 91/576+7/192*sqrt(21)-(585/1568*b[8])*sqrt(21)-405/224*b[8]; A[ 9, 8] := (585/1568*b[8])*sqrt(21)+405/224*b[8]; A[10, 1] := 1/14; A[10, 2] := 0; A[10, 3] := -(6/49*a[10, 5])*sqrt(21)-2/7*a[10,5]+2/147*sqrt(21)+2/63; A[10, 4] := 1/9-a[10, 5]; A[10, 5] := a[10, 5]; A[10, 6] := 2/7*a[10, 5]-803/2205+(6/49*a[10,5])*sqrt(21)-59/735*sqrt(21); A[10, 7] := 1/9+(1/42)*sqrt(21)+2295/686*b[8]+(495/686*b[8])*sqrt(21); A[10, 8] := -2295/686*b[8]-(495/686*b[8])*sqrt(21); A[10, 9] := 4/35*21^(1/2)+132/245; A[11, 1] := 0; A[11, 2] := 0; A[11, 3] := (2/3*a[10, 5])*sqrt(21)-2/3*a[10, 5]-2/27*sqrt(21)+2/27; A[11, 4] := -(7/6*a[10, 5])*sqrt(21)+7/54*sqrt(21)+49/6*a[10,5]-49/54; A[11, 5] := 7/27*sqrt(21)-77/54-49/6*a[10, 5]+(7/6*a[10,5])*sqrt(21); A[11, 6] := 2/3*a[10, 5]-64/135-(2/3*a[10,5])*sqrt(21)+94/135*sqrt(21); A[11, 7] := 7/18-(265/98*b[8])*sqrt(21)-215/14*b[8]; A[11, 8] := (265/98*b[8])*sqrt(21)+215/14*b[8]; A[11, 9] := -28/45*sqrt(21)+28/45; A[11, 10] := 49/18-7/18*sqrt(21); # b[i]: A[12, 1] := 1/20; A[12, 2] := 0; A[12, 3] := 0; A[12, 4] := 0; A[12, 5] := 0; A[12, 6] := 0; A[12, 7] := -b[8]+49/180; A[12, 8] := b[8]; A[12, 9] := 16/45; A[12, 10] := 49/180; A[12, 11] := 1/20; simplify(A); end proc: