Sponsored links: Algebra eBooks
 

Related

error-load-matrix

load(lsquares);

D: matrix( [3.2, 12,...

ti:22;

Calculate

error-exp-sqrt

lambda_Z:3;

lambda_X:12;

N:50;

Calculate

error-exp-sqrt

lambda_Z:3;

lambda_X:12;

N:50;

Calculate

error-exp-sqrt

lambda_Z:3;

lambda_X:12;

N:50;

Calculate

error-exp-integrate-load-sqrt

load(dblint) ;

lambda_Z:3 ;

lambda_X:12 ;

Calculate

error-exp-sqrt

lambda_Z:3;

lambda_X:12;

N:50;

Calculate

error-exp-sqrt

lambda_Z:3;

lambda_X:12;

N:50;

Calculate

error-exp-sqrt

lambda_Z:3;

lambda_X:12;

N:50;

Calculate

error-load-matrix

load(lsquares);

D: matrix( [3.2, 12,...

ti:22;

Calculate

error

Run Example
(%i1)load(lsquares);
(%o1)         /usr/share/maxima/5.21.1/share/contrib/lsquares.mac
(%i2) D: matrix(  [3.2, 12, 60],  [3.2, 4, 89],  [3.2, -4, 115],  [3.2, -8, 126],  [3.2, 0, 102],  [3.2, 8, 75],  [3.2, 12, 60],  [2, -16, 99],  [2, -12, 93],  [2, -8, 86],  [2, -4, 78.5],  [2, 0, 71],  [2, 4, 63],  [2, 8, 55],  [2, 12, 46],  [1, -16, 59],  [1, -12, 56],  [1, -8, 54],  [1, -4, 50],  [1, 0, 46],  [1, 4, 42],  [1, 8, 39],  [1, 12, 34]);
                              [ 3.2   12    60  ]
                              [                 ]
                              [ 3.2   4     89  ]
                              [                 ]
                              [ 3.2  - 4   115  ]
                              [                 ]
                              [ 3.2  - 8   126  ]
                              [                 ]
                              [ 3.2   0    102  ]
                              [                 ]
                              [ 3.2   8     75  ]
                              [                 ]
                              [ 3.2   12    60  ]
                              [                 ]
                              [  2   - 16   99  ]
                              [                 ]
                              [  2   - 12   93  ]
                              [                 ]
                              [  2   - 8    86  ]
                              [                 ]
                              [  2   - 4   78.5 ]
                              [                 ]
(%o2)                         [  2    0     71  ]
                              [                 ]
                              [  2    4     63  ]
                              [                 ]
                              [  2    8     55  ]
                              [                 ]
                              [  2    12    46  ]
                              [                 ]
                              [  1   - 16   59  ]
                              [                 ]
                              [  1   - 12   56  ]
                              [                 ]
                              [  1   - 8    54  ]
                              [                 ]
                              [  1   - 4    50  ]
                              [                 ]
                              [  1    0     46  ]
                              [                 ]
                              [  1    4     42  ]
                              [                 ]
                              [  1    8     39  ]
                              [                 ]
                              [  1    12    34  ]
(%i3) E:matrix(  [0.4, -16, 42],  [0.4, -12, 41],  [0.4, -8, 39],  [0.4, -4, 37],  [0.4, 0, 35],  [0.4, 4, 33],  [0.4, 8, 31],  [0.4, 12, 29]);
                               [ 0.4  - 16  42 ]
                               [               ]
                               [ 0.4  - 12  41 ]
                               [               ]
                               [ 0.4  - 8   39 ]
                               [               ]
                               [ 0.4  - 4   37 ]
(%o3)                          [               ]
                               [ 0.4   0    35 ]
                               [               ]
                               [ 0.4   4    33 ]
                               [               ]
                               [ 0.4   8    31 ]
                               [               ]
                               [ 0.4   12   29 ]
(%i4) D : addrow(D, E);
                              [ 3.2   12    60  ]
                              [                 ]
                              [ 3.2   4     89  ]
                              [                 ]
                              [ 3.2  - 4   115  ]
                              [                 ]
                              [ 3.2  - 8   126  ]
                              [                 ]
                              [ 3.2   0    102  ]
                              [                 ]
                              [ 3.2   8     75  ]
                              [                 ]
                              [ 3.2   12    60  ]
                              [                 ]
                              [  2   - 16   99  ]
                              [                 ]
                              [  2   - 12   93  ]
                              [                 ]
                              [  2   - 8    86  ]
                              [                 ]
                              [  2   - 4   78.5 ]
                              [                 ]
                              [  2    0     71  ]
                              [                 ]
                              [  2    4     63  ]
                              [                 ]
                              [  2    8     55  ]
                              [                 ]
                              [  2    12    46  ]
                              [                 ]
(%o4)                         [  1   - 16   59  ]
                              [                 ]
                              [  1   - 12   56  ]
                              [                 ]
                              [  1   - 8    54  ]
                              [                 ]
                              [  1   - 4    50  ]
                              [                 ]
                              [  1    0     46  ]
                              [                 ]
                              [  1    4     42  ]
                              [                 ]
                              [  1    8     39  ]
                              [                 ]
                              [  1    12    34  ]
                              [                 ]
                              [ 0.4  - 16   42  ]
                              [                 ]
                              [ 0.4  - 12   41  ]
                              [                 ]
                              [ 0.4  - 8    39  ]
                              [                 ]
                              [ 0.4  - 4    37  ]
                              [                 ]
                              [ 0.4   0     35  ]
                              [                 ]
                              [ 0.4   4     33  ]
                              [                 ]
                              [ 0.4   8     31  ]
                              [                 ]
                              [ 0.4   12    29  ]
(%i5) ti:22;
(%o5)                                 22
(%i6) tamin:-30;
(%o6)                                - 30
(%i7) r:lsquares_estimates(D, [m, ta, tvl], tvl = ti+(m*b^c)*(((ti-ta)/(ti-tamin))^(a)), [a, b,c]);
*************************************************
  N=    3   NUMBER OF CORRECTIONS=25
       INITIAL VALUES
 F=  2.101931255964879D+03   GNORM=  9.218838988104815D+01
*************************************************

   I  NFN     FUNC                    GNORM                   STEPLENGTH

   1    7     1.116228467499764D+03   5.181482092741437D+02   3.792158188175497D-02  
   2   11     1.115139217279903D+03   4.517933553619771D+02   6.156039739677396D-03  
   3   15     3.772201212301974D+02   1.464357581324870D+03   1.279783986132719D+01  
   4   17     7.783971751722766D+01   1.150713443659835D+03   2.039382076736760D-01  
   5   19     5.003402304305900D+00   5.919845697246858D+01   1.995821879972414D-01  
   6   20     4.211639387760370D+00   5.820407389162347D+01   1.000000000000000D+00  
   7   21     4.088456788513105D+00   5.277714141954149D+00   1.000000000000000D+00  
   8   22     4.087286836630353D+00   2.094011852767751D-02   1.000000000000000D+00  
   9   23     4.087286523816272D+00   3.178623521461571D-05   1.000000000000000D+00  

 THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.
 IFLAG = 0
(%o7) [[a = 0.87205051345354, b = 3.011617383235275, c = 3.585922792420821]]
(%i8) a: rhs(r[1][1]);
(%o8)                          0.87205051345354
(%i9) b: rhs(r[1][2]);
(%o9)                          3.011617383235275
(%i10) c: rhs(r[1][3]);
(%o10)                         3.585922792420821
(%i11) error(D) := (ti+(col(D, 1)*b)*((ti-col(D, 2))/(ti-tamin))^^(1/a+c*col(D, 1))) - col(D,3);
                                       ti - col(D, 2) <1/a + c col(D, 1)>
(%o11) error(D) := ti + (col(D, 1) b) (--------------)
                                         ti - tamin
                                                                    - col(D, 3)
(%i12) e(m, ta, tvl) := ti+(m*b^c)*((ti-ta)/(ti-tamin))^(a) - tvl;
                                         c    ti - ta   a
(%o12)         e(m, ta, tvl) := ti + (m b ) (----------)  - tvl
                                             ti - tamin
(%i13) for i:1 while i <
= length(first(transpose(D))) do disp (  e(D[i][1], D[i][2], D[i][3]));
                               1.600256016752056

                              - 0.88370438260647

                              - 1.887971702076186

                              - 0.77810475937126

                              - 1.239608130221058

                               0.10421852552332

                               1.600256016752056

                               2.282729637294366

                               0.95392240560348

                               0.51368452539296

                               0.44501768620238

                               0.22524491861184

                               0.32268476087096

                               0.19013657845207

                               0.75016001047003

                               2.641364818647183

                               1.976961202801739

                               0.25684226269648

                               0.47250884310119

                               0.61262245930592

                               0.66134238043548

                              - 0.40493171077397

                               0.37508000523501

                              - 4.143454072541125

                              - 4.609215518879303

                              - 4.097263094921409

                              - 3.610996462759523

                              - 3.154951016277632

                              - 2.73546304782581

                              - 2.361972684309585

                              - 2.049967997905995

(%o13)                               done
(%i14) 
Run Example
get_deg_occ(eq,alist):=block(	[atomlistlen,deglist,occlist,i],	atomlistlen: length(alist),	deglist: makelist(0,ii,1,atomlistlen),	occlist: makelist(0,ii,1,atomlistlen),	for i:1 thru atomlistlen do (		deglist[i]: hipow(eq,alist[i]),		if deglist[i]>
0 then (			occlist[i]: 1		)	),	return([deglist,occlist]));
(%o1) get_deg_occ(eq, alist) := block([atomlistlen, deglist, occlist, i], 
atomlistlen : length(alist), deglist : makelist(0, ii, 1, atomlistlen), 
occlist : makelist(0, ii, 1, atomlistlen), 
for i thru atomlistlen do (deglist  : hipow(eq, alist ), 
                                  i                  i
if deglist  > 0 then occlist  : 1), return([deglist, occlist]))
          i                 i
(%i2) remove_col_row(dmatrix,omatrix,aelist,atomlistlen,eqlistlen,sumvalue):=block(	[degmatrix,occmatrix,llist,llistlen,latomlistlen,leqlistlen,ret,pos,i],	degmatrix: copymatrix(dmatrix),	occmatrix: copymatrix(omatrix),	llist: copylist(aelist),	latomlistlen: atomlistlen,	leqlistlen: eqlistlen,	ret: false,	while ret=false do (		pos: 0,		for i:1 thru latomlistlen do (			if sum(occmatrix[ii,i],ii,1,leqlistlen)=sumvalue then (				if sumvalue=0 then (					print("Column ",i," has sum=0 (atom=",llist[i],") and must be removed")				)				else (					/* sumvalue=1 */					print("Column ",i," has sum=1 and must be removed")				),				pos: i,				return(pos)			)		),		if pos>
0 then (			if sumvalue=0 then (				/* remove column pos from degmatrix and occmatrix, remove pos-th atom from atomlist */				degmatrix: submatrix(degmatrix,pos),				occmatrix: submatrix(occmatrix,pos),				llist: delete(llist[pos],llist),				latomlistlen: length(llist)			)			else (				/* sumvalue=1 */				/* remove row pos from degmatrix and occmatrix, remove pos-th equation from eqlist */				degmatrix: submatrix(pos,degmatrix),				occmatrix: submatrix(pos,occmatrix),				llist: delete(llist[pos],llist),				leqlistlen: length(llist)			)		)		else (			ret: true		),		if sumvalue=0 and latomlistlen=0 then (			error("no atoms remaining")		)		else if sumvalue=1 and leqlistlen=0 then (			error("no equations remaining")		)	),	return([degmatrix,occmatrix,llist]));
(%o2) remove_col_row(dmatrix, omatrix, aelist, atomlistlen, eqlistlen, 
sumvalue) := block([degmatrix, occmatrix, llist, llistlen, latomlistlen, 
leqlistlen, ret, pos, i], degmatrix : copymatrix(dmatrix), 
occmatrix : copymatrix(omatrix), llist : copylist(aelist), 
latomlistlen : atomlistlen, leqlistlen : eqlistlen, ret : false, 
while ret = false do (pos : 0, for i thru latomlistlen 
do if sum(occmatrix     , ii, 1, leqlistlen) = sumvalue
                   ii, i
 then (if sumvalue = 0 then print("Column ", i, " has sum=0 (atom=", llist , 
                                                                          i
") and must be removed") else print("Column ", i, 
" has sum=1 and must be removed"), pos : i, return(pos)), 
if pos > 0 then (if sumvalue = 0 then (degmatrix : submatrix(degmatrix, pos), 
occmatrix : submatrix(occmatrix, pos), llist : delete(llist   , llist), 
                                                           pos
latomlistlen : length(llist)) else (degmatrix : submatrix(pos, degmatrix), 
occmatrix : submatrix(pos, occmatrix), llist : delete(llist   , llist), 
                                                           pos
leqlistlen : length(llist))) else ret : true, 
if (sumvalue = 0) and (latomlistlen = 0) then error("no atoms remaining")
 else (if (sumvalue = 1) and (leqlistlen = 0)
 then error("no equations remaining"))), return([degmatrix, occmatrix, llist]))
(%i3) pdesys2pde(alist,atlist,elist,etlist):=block(	[atomlist,atomtlist,eqlist,eqtlist,atomlistlen,atomtlistlen,eqlistlen,eqtlistlen,degmatrix,occmatrix,deglist,occlist,i,j,deg,ret,pos,neqlist,neqlistlen,nneqlist],	eqlist: copylist(elist),	eqtlist: copy(etlist),	atomlist: copylist(alist),	atomtlist: copylist(atlist),	eqlistlen: length(eqlist),	eqtlistlen: length(eqtlist),	atomlistlen: length(atomlist),	atomtlistlen: length(atomtlist),	/* for each equation, get the highest degree for each t-atom */	degmatrix: matrix(),	/* for each atom, get the number of occurrences in the equation set */	occmatrix: matrix(),	for i:1 thru eqtlistlen do (		[deglist,occlist]: get_deg_occ(eqtlist[i],atomtlist),		degmatrix: addrow(degmatrix,deglist),		occmatrix: addrow(occmatrix,occlist)	),print("Start: ",degmatrix),print("Start: ",occmatrix),return(0),	/* remove (recursively) all rows, whose sum is 1, that is, all equation containing atom occuring only once */	[degmatrix,occmatrix,eqlist]: remove_col_row(degmatrix,occmatrix,eqlist,atomlistlen,eqlistlen,1),	eqlistlen: length(eqlist),/*	ret: false,	while ret=false do (		pos: 0,		for i:1 thru atomlistlen do (			if sum(occmatrix[ii,i],ii,1,eqlistlen)=1 then (				print("Column ",i," has sum=1 and must be removed"),				pos: i,				return(pos)			)		),		if pos>
0 then (			-- remove row pos from degmatrix and occmatrix, remove pos-th equation from eqlist --			degmatrix: submatrix(pos,degmatrix),			occmatrix: submatrix(pos,occmatrix),			eqlist: delete(eqlist[pos],eqlist),			eqlistlen: length(eqlist)		)		else (			ret: true		),		if eqlistlen=0 then (			error("no equations remaining")		)	),*/	/* now remove (recursively) all columns, whose sum is 0: these variables are not used */	[degmatrix,occmatrix,atomlist]: remove_col_row(degmatrix,occmatrix,atomlist,atomlistlen,eqlistlen,0),	atomlistlen: length(atomlist),/*	ret: false,	while ret=false do (		pos: 0,		for i:1 thru atomlistlen do (			if sum(occmatrix[ii,i],ii,1,eqlistlen)=0 then (				print("Column ",i," has sum=0 (atom=",atomlist[i],") and must be removed"),				pos: i,				return(pos)			)		),		if pos>
0 then (			/* remove column pos from degmatrix and occmatrix, remove pos-th atom from atomlist */			degmatrix: submatrix(degmatrix,pos),			occmatrix: submatrix(occmatrix,pos),			atomlist: delete(atomlist[pos],atomlist),			atomlistlen: length(atomlist)		)		else (			ret: true		),		if atomlistlen=0 then (			error("no atoms remaining")		)	),*//*	for i:atomlistlen thru 1 step -1 do (		if sum(occmatrix[ii,i],ii,1,eqlistlen)=0 then (			-- remove column i from degmatrix and occmatrix, remove i-th atom from atomlist --			degmatrix: submatrix(degmatrix,i),			occmatrix: submatrix(occmatrix,i),			print("Column ",i," has sum=0 (atom=",atomlist[i],") and must be removed"),			atomlist: delete(atomlist[pos],atomlist)		)	),	atomlistlen: length(atomlist),*/print("End: ",atomlist),print("End: ",eqlist),print("End: ",degmatrix),print("End: ",occmatrix),return(0),	/* now, each atom exists in at least 2 equations */	/* reduce the set of equations, by eliminating the atoms starting by the highest derivative (x+t) */print("***"),	/* WWWHHH */	pos: atomlistlen,	neqlist: [],	for i:eqlistlen thru 1 step -1 do (		if occmatrix[i,pos]=1 then (print(i),			neqlist: append(neqlist,[eqlist[i]]),			/* remove row i from degmatrix and occmatrix */			degmatrix: submatrix(degmatrix,i),			occmatrix: submatrix(occmatrix,i)		)	),	nneqlist: [],	/* choose the data based on the degree! */	/* this is the basic variation: first equation is fixed */	for i:2 thru neqlistlen do (		nneqlist: append(nneqlist,[resultant(neqlist[1],neqlist[i],atomlist[pos])])	),	neqlist: copylist(nneqlist),	neqlistlen: length(neqlist),	atomlist: delete(atomlist[pos],atomlist),	atomlistlen: length(atomlist),	/* add row to degmatrix and occmatrix */	for i:1 thru neqlistlen do (		[deglist,occlist]: get_deg_occ(neqlist[i],atomlist),		degmatrix: addrow(degmatrix,deglist),		occmatrix: addrow(occmatrix,occlist)	),	/* now remove columns with sum=0 */	/* WWWHHH */	return(0));
(%o3) pdesys2pde(alist, atlist, elist, etlist) := 
block([atomlist, atomtlist, eqlist, eqtlist, atomlistlen, atomtlistlen, 
eqlistlen, eqtlistlen, degmatrix, occmatrix, deglist, occlist, i, j, deg, ret, 
pos, neqlist, neqlistlen, nneqlist], eqlist : copylist(elist), 
eqtlist : copy(etlist), atomlist : copylist(alist), 
atomtlist : copylist(atlist), eqlistlen : length(eqlist), 
eqtlistlen : length(eqtlist), atomlistlen : length(atomlist), 
atomtlistlen : length(atomtlist), degmatrix : matrix(), occmatrix : matrix(), 
for i thru eqtlistlen do ([deglist, occlist] : 
get_deg_occ(eqtlist , atomtlist), degmatrix : addrow(degmatrix, deglist), 
                   i
occmatrix : addrow(occmatrix, occlist)), print("Start: ", degmatrix), 
print("Start: ", occmatrix), return(0), 
[degmatrix, occmatrix, eqlist] : remove_col_row(degmatrix, occmatrix, eqlist, 
atomlistlen, eqlistlen, 1), eqlistlen : length(eqlist), 
[degmatrix, occmatrix, atomlist] : remove_col_row(degmatrix, occmatrix, 
atomlist, atomlistlen, eqlistlen, 0), atomlistlen : length(atomlist), 
print("End: ", atomlist), print("End: ", eqlist), print("End: ", degmatrix), 
print("End: ", occmatrix), return(0), print("***"), pos : atomlistlen, 
neqlist : [], for i from eqlistlen step - 1 thru 1 
do if occmatrix       = 1 then (print(i), 
               i, pos
neqlist : append(neqlist, [eqlist ]), degmatrix : submatrix(degmatrix, i), 
                                 i
occmatrix : submatrix(occmatrix, i)), nneqlist : [], 
for i from 2 thru neqlistlen do nneqlist : 
append(nneqlist, [resultant(neqlist , neqlist , atomlist   )]), 
                                   1         i          pos
neqlist : copylist(nneqlist), neqlistlen : length(neqlist), 
atomlist : delete(atomlist   , atomlist), atomlistlen : length(atomlist), 
                          pos
for i thru neqlistlen do ([deglist, occlist] : 
get_deg_occ(neqlist , atomlist), degmatrix : addrow(degmatrix, deglist), 
                   i
occmatrix : addrow(occmatrix, occlist)), return(0))
(%i4) pdesys2pde([v,vx,vxx,vxxx],[vt,vxt],[vt-ux*vxx,utt-vt*vx-v*vxt,vxt-uxx*vxx-ux*vxxx],[ut-v*vx,uxt-vx^2-v*vxx,uxxt-3*vx*vxx-v*vxxx]);
        [ 0  0 ]
        [      ]
Start:  [ 0  0 ] 
        [      ]
        [ 0  0 ]
        [ 0  0 ]
        [      ]
Start:  [ 0  0 ] 
        [      ]
        [ 0  0 ]
(%o4)                                  0
(%i5) 
Run Example
sphere_error:(x[i]-a)^2 + (y[i]-b)^2 + (z[i]-c)^2 - r2;
                                  2           2           2
(%o1)              - r2 + (z  - c)  + (y  - b)  + (x  - a)
                            i           i           i
(%i2) sphere_error_square:sphere_error^2;
                                  2           2           2 2
(%o2)             (- r2 + (z  - c)  + (y  - b)  + (x  - a) )
                            i           i           i
(%i3) declare(sum,linear);
(%o3)                                done
(%i4) sphere_error_sum:sum(sphere_error_square, i, 1, n);
                n
               ====
               \                     2           2           2 2
(%o4)           >    (- r2 + (z  - c)  + (y  - b)  + (x  - a) )
               /               i           i           i
               ====
               i = 1
(%i5) sphere_da: diff(sphere_error_sum, a);
              n
             ====
             \                              2           2           2
(%o5)    - 4  >    (x  - a) (- r2 + (z  - c)  + (y  - b)  + (x  - a) )
             /       i                i           i           i
             ====
             i = 1
(%i6) sphere_db: diff(sphere_error_sum, b);
              n
             ====
             \                              2           2           2
(%o6)    - 4  >    (y  - b) (- r2 + (z  - c)  + (y  - b)  + (x  - a) )
             /       i                i           i           i
             ====
             i = 1
(%i7) sphere_dc: diff(sphere_error_sum, c);
              n
             ====
             \                              2           2           2
(%o7)    - 4  >    (z  - c) (- r2 + (z  - c)  + (y  - b)  + (x  - a) )
             /       i                i           i           i
             ====
             i = 1
(%i8) sphere_dr2: diff(sphere_error_sum, r2);
                     n                 n                 n
                    ====              ====              ====
                    \             2   \             2   \             2
(%o8) - 2 (- n r2 +  >    (z  - c)  +  >    (y  - b)  +  >    (x  - a) )
                    /       i         /       i         /       i
                    ====              ====              ====
                    i = 1             i = 1             i = 1
(%i9) sphere_da_exp: expand(sphere_da=0);
                       n
                      ====
                      \                   2          2        3
(%o9) - 4 a n r2 + 4 ( >    x ) r2 + 4 a c  n + 4 a b  n + 4 a  n
                      /      i
                      ====
                      i = 1
      n                 n              n                   n
     ====              ====           ====                ====
     \         2       \      2       \                   \
 - 4  >    x  z  + 4 a  >    z  + 8 c  >    x  z  - 8 a c  >    z
     /      i  i       /      i       /      i  i         /      i
     ====              ====           ====                ====
     i = 1             i = 1          i = 1               i = 1
      n                 n              n                   n            n
     ====              ====           ====                ====         ====
     \         2       \      2       \                   \            \      3
 - 4  >    x  y  + 4 a  >    y  + 8 b  >    x  y  - 8 a b  >    y  - 4  >    x
     /      i  i       /      i       /      i  i         /      i     /      i
     ====              ====           ====                ====         ====
     i = 1             i = 1          i = 1               i = 1        i = 1
         n               n               n                n
        ====            ====            ====             ====
        \      2      2 \             2 \              2 \
 + 12 a  >    x  - 4 c   >    x  - 4 b   >    x  - 12 a   >    x  = 0
        /      i        /      i        /      i         /      i
        ====            ====            ====             ====
        i = 1           i = 1           i = 1            i = 1
(%i10) sphere_db_exp: expand(sphere_db=0);
                        n
                       ====
                       \                   2        3        2
(%o10) - 4 b n r2 + 4 ( >    y ) r2 + 4 b c  n + 4 b  n + 4 a  b n
                       /      i
                       ====
                       i = 1
      n                 n              n                   n            n
     ====              ====           ====                ====         ====
     \         2       \      2       \                   \            \      3
 - 4  >    y  z  + 4 b  >    z  + 8 c  >    y  z  - 8 b c  >    z  - 4  >    y
     /      i  i       /      i       /      i  i         /      i     /      i
     ====              ====           ====                ====         ====
     i = 1             i = 1          i = 1               i = 1        i = 1
         n            n                 n                  n
        ====         ====              ====               ====
        \      2     \      2          \                2 \
 + 12 b  >    y  - 4  >    x  y  + 8 a  >    x  y  - 4 c   >    y
        /      i     /      i  i       /      i  i        /      i
        ====         ====              ====               ====
        i = 1        i = 1             i = 1              i = 1
          n               n              n                n
         ====            ====           ====             ====
       2 \             2 \              \      2         \
 - 12 b   >    y  - 4 a   >    y  + 4 b  >    x  - 8 a b  >    x  = 0
         /      i        /      i       /      i         /      i
         ====            ====           ====             ====
         i = 1           i = 1          i = 1            i = 1
(%i11) sphere_dc_exp: expand(sphere_dc=0);
                        n                                               n
                       ====                                            ====
                       \                 3        2          2         \      3
(%o11) - 4 c n r2 + 4 ( >    z ) r2 + 4 c  n + 4 b  c n + 4 a  c n - 4  >    z
                       /      i                                        /      i
                       ====                                            ====
                       i = 1                                           i = 1
         n            n                 n               n
        ====         ====              ====            ====
        \      2     \      2          \               \      2
 + 12 c  >    z  - 4  >    y  z  + 8 b  >    y  z  - 4  >    x  z
        /      i     /      i  i       /      i  i     /      i  i
        ====         ====              ====            ====
        i = 1        i = 1             i = 1           i = 1
        n                   n               n               n
       ====                ====            ====            ====
       \                 2 \             2 \             2 \
 + 8 a  >    x  z  - 12 c   >    z  - 4 b   >    z  - 4 a   >    z
       /      i  i         /      i        /      i        /      i
       ====                ====            ====            ====
       i = 1               i = 1           i = 1           i = 1
        n                n              n                n
       ====             ====           ====             ====
       \      2         \              \      2         \
 + 4 c  >    y  - 8 b c  >    y  + 4 c  >    x  - 8 a c  >    x  = 0
       /      i         /      i       /      i         /      i
       ====             ====           ====             ====
       i = 1            i = 1          i = 1            i = 1
(%i12) augcoefmatrix([sphere_da_subst, sphere_db_subst, sphere_dc_subst], [a, b, c]);
                         [ 0  0  0  sphere_da_subst ]
                         [                          ]
(%o12)                   [ 0  0  0  sphere_db_subst ]
                         [                          ]
                         [ 0  0  0  sphere_dc_subst ]
(%i13) 

Related Help

Help for Error