| blob | 3a8915c41ca260d7b1527547ddfac6893df0a021 |
1 /*************** -*- Mode: MACSYMA; Package: MAXIMA -*- ******************/
3 (kill(all), done);
4 done;
6 /* Examples from Beyer (1984), an early paper on algsys, which references
7 an internal General Motors report by Morgan (1981). Morgan (1983) has
8 the same title and contains some of the problems.
10 Beyer, William A., Solution of Simultaneous Polynomial Equations by
11 Elimination in MACSYMA, Proceedings of the 1984 MACSYMA Users' Conference,
12 p 110-120
14 Morgan, Alexander P., A Method for Computing All Solutions to Systems of
15 Polynomials Equations, ACM Transactions on Mathematical Software,
16 9:1-17 (1983)
17 */
19 /* Beyer example 1, p117 */
20 algsys([x^2+x*y-1,y^2+x-5],[x,y]);
21 [[x = -3.174087266251113,y = 2.859036144578313],
22 [x = 2.152189020152884,y = -1.687545787545788],
23 [x = 0.3937104979189148,y = 2.146226811046756],
24 [x = -0.3718122830843519,y = -2.317717206132879]];
26 /* Beyer example 2, p117, Morgan problem 1 */
27 algsys([4*x^3-3*x-y,x^2-y],[x,y]);
28 [[x = 1,y = 1],[x = -3/4,y = 9/16],[x = 0,y = 0]];
30 /* Beyer example 3, p117, Morgan problem 2 */
31 algsys([4*(x+y),4*(x+y)+(x-y)*((x-2)^2+y^2-1)],[x,y]);
32 [[x = -((sqrt(2)*%i-2)/2), y = (sqrt(2)*%i-2)/2],
33 [x = (sqrt(2)*%i+2)/2, y = -((sqrt(2)*%i+2)/2)],
34 [x = 0,y = 0]];
36 /* Beyer example 4, p117. Don't understand notation in paper */
38 /* Beyer example 5, p118. 2 real and 6 complex solutions */
39 algsys([
40 x^2+2*y^2-4,
41 x^2+y^2+z-8,
42 (x-1)^2+(2*y-sqrt(2))^2+(z-5)^2-4],
43 [x,y,z]);
44 [[x = 0, y = sqrt(2), z = 6],
45 [x = 2.998824877038494 - 0.879551236563151*%i,
46 y = (- 1.656382492469821*%i) - 0.796198988104073,
47 z = 2.637620128835544*%i + 1.890329867297412],
48 [x = 0.879551236563151*%i + 2.998824877038494,
49 y = 1.656382492469821*%i - 0.7961989881040729,
50 z = 1.890329867297412 - 2.637620128835543*%i],
51 [x = (- 1.511881394988993*%i) - 2.77266624385368,
52 y = 1.573386231515891*%i - 1.332140330399878,
53 z = 3.299053626354102 - 4.191942508596404*%i],
54 [x = 1.511881394988993*%i - 2.77266624385368,
55 y = (- 1.573386231515895*%i) - 1.332140330399874,
56 z = 4.191942508596393*%i + 3.299053626354096],
57 [x = (- 1.06051780109088*%i) - 1.226158633184815,
58 y = 1.421232537317398 - 0.4574772330741239*%i,
59 z = 5.81061650634849 - 1.30036305745376*%i],
60 [x = 1.06051780109088*%i - 1.226158633184815,
61 y = 0.4574772330741241*%i + 1.421232537317397,
62 z = 1.300363057453758*%i + 5.810616506348491],
63 [x = 2, y = 0, z = 4]];
65 /* Beyer example 6, p118 */
66 block([%rnum:0],
67 algsys([x+10*y,z+w,(z-2*z)^2,(z-w)^2],[w,x,y,z]));
68 [[w = 0, x = %r1 , y = -(%r1/10), z = 0]];
70 /* Beyer example 7, p118. Supposed to be two real solutions.
71 This solution seems correct. Typo in paper? */
72 algsys([
73 x+y+z+w-1,
74 x+y-z+w-3,
75 x^2+y^2+z^2+w^2-4,
76 (x-1)^2+y^2+z^2+w^2],
77 [w,x,y,z]);
78 [[w = -((3^(3/2)*%i+1)/4), x = 5/2, y = (3^(3/2)*%i-1)/4, z = -1],
79 [w = (3^(3/2)*%i-1)/4, x = 5/2, y = -((3^(3/2)*%i+1)/4), z = -1]];
81 /* Check Beyer example 7. Preprocess with poly_reduced_grobner(). OK */
82 algsys([z+1,2*x-5,4*y^2+2*y+7,2*y+2*w+1],[w,x,y,z]);
83 [[w = -((3^(3/2)*%i+1)/4), x = 5/2, y = (3^(3/2)*%i-1)/4, z = -1],
84 [w = (3^(3/2)*%i-1)/4, x = 5/2, y = -((3^(3/2)*%i+1)/4), z = -1]];
86 /* Beyer example 8, p118. Two real solutions */
87 algsys([
88 x1^2+x2^2+x3^2+x4-6,
89 x1+x2+x5+x6-1,
90 x1+3*x3+x6-3,
91 x1+x2+x4-2,
92 2*x3+x6,
93 x5+x6],
94 [x1,x2,x3,x4,x5,x6]);
95 [[x1 = 1, x2 = 0, x3 = 2, x4 = 1, x5 = 4, x6 = -4],
96 [x1 = 5/3,x2 = -2/3,x3 = 4/3,x4 = 1, x5 = 8/3, x6 = -8/3]];
98 /* Beyer example 9, p118. Morgan problem 4
99 Three real and two complex solutions */
100 algsys(
101 [x1+x1+x2+x3+x4+x5-6=0,
102 x2+x1+x2+x3+x4+x5-6=0,
103 x3+x1+x2+x3+x4+x5-6=0,
104 x4+x1+x2+x3+x4+x5-6=0,
105 x1*x2*x3*x4*x5-1=0],
106 [x1,x2,x3,x4,x5]);
107 [[x1 = 0.916354556803995,
108 x2 = 0.916354556803995,
109 x3 = 0.916354556803995,
110 x4 = 0.916354556803995,
111 x5 = 1.418227215980025],
112 [x1 = -0.5790430889759374,
113 x2 = -0.5790430889759374,
114 x3 = -0.5790430889759374,
115 x4 = -0.5790430889759374,
116 x5 = 8.895215311004785],
117 [x1 = (-0.6100917318163317*%i)-0.06865574701986676,
118 x2 = (-0.6100917318163317*%i)-0.06865574701986676,
119 x3 = (-0.6100917318163317*%i)-0.06865574701986676,
120 x4 = (-0.6100917318163317*%i)-0.06865574701986676,
121 x5 = 3.050458659081659*%i+6.343278735099334],
122 [x1 = 0.6100917318163317*%i-0.06865574701986676,
123 x2 = 0.6100917318163318*%i-0.06865574701986678,
124 x3 = 0.6100917318163317*%i-0.06865574701986676,
125 x4 = 0.6100917318163317*%i-0.06865574701986676,
126 x5 = 6.343278735099333-3.050458659081658*%i],
127 [x1 = 1, x2 = 1, x3 = 1, x4 = 1, x5 = 1]];
129 /* Beyer example 10, p118. 4 real one-parameter solutions */
130 block([%rnum:0],
131 algsys([x^2+y^2-1,x^2+y^2+z^2-5=0],[x,y,z]));
132 [[x = %r1, y = -sqrt(1-%r1^2), z = 2],
133 [x = %r2, y = sqrt(1-%r2^2), z = 2],
134 [x = %r3, y = -sqrt(1-%r3^2), z = -2],
135 [x = %r4, y = sqrt(1-%r4^2), z = -2]];
137 /* Morgan problem 3. Single solution (0,0,0,0) */
138 algsys([x+10*y,sqrt(5)*(z-w),(y-2*z)^2,sqrt(10)*(x-w)^2],[w,x,y,z]);
139 [[w = 0, x = 0, y = 0, z = 0]];
141 /* SF bug #3208 Can't solve this simple equation */
142 algsys([x^2+y^2=1,(x-2)^2+(y-3)^2=9],[x,y]);
143 [[x = -((3^(5/2)-10)/26),y = (2*3^(3/2)+15)/26],
144 [x = (3^(5/2)+10)/26,y = -((2*3^(3/2)-15)/26)]];
146 /* Reduced from SF bug #3208 above */
147 algsys([676*y^2-20*3^(5/2)-333,676*y^2-4056*y+28*3^(7/2)+2007],[y]);
148 [[y = (2*3^(3/2)+15)/26]];
150 /* Reduced from SF bug #3208 above */
151 algsys([676*y^2+20*3^(5/2)-333,676*y^2-4056*y-28*3^(7/2)+2007],[y]);
152 [[y = -((2*3^(3/2)-15)/26)]];
154 /* from example(solve) */
155 algsys([4*x^2-y^2=12,x*y-x=2],[x,y]);
156 [[x = 2, y = 2],
157 [x = 0.5202594388652008*%i - 0.1331240357358706,
158 y = 0.07678378523787788 - 3.608003221870287*%i],
159 [x = -0.5202594388652008*%i - 0.1331240357358706,
160 y = 3.608003221870287*%i + 0.07678378523787788],
161 [x = -1.733751846381093, y = -0.1535675710019696]];
163 /* failures from example(algsys) during testing */
164 block( [
165 %rnum:0,
166 f1:2*x*(1-l1)-2*(x-1)*l2,
167 f2:l2-l1,
168 f3:l1*(1-x**2-y),
169 f4:l2*(y-(x-1)**2)
170 ],
171 algsys([f1,f2,f3,f4],[x,y,l1,l2]));
172 [[x = 0, y = %r1, l1 = 0, l2 = 0],
173 [x = 1, y = 0, l1 = 1, l2 = 1]];
175 block( [
176 f1:x**2-y**2,
177 f2:x**2-x+2*y**2-y-1
178 ],
179 algsys([f1,f2],[x,y]));
180 [[x = -(1/sqrt(3)), y = 1/sqrt(3)],
181 [x = 1/sqrt(3), y = -(1/sqrt(3))],
182 [x = -(1/3),y = -(1/3)],
183 [x = 1,y = 1]];
185 /* failure in example(charpoly) during testing */
186 algsys([x2-2*x1=0,x2^2+x1^2=1],[x1,x2]);
187 [[x1 = -(1/sqrt(5)), x2 = -(2/sqrt(5))],
188 [x1 = 1/sqrt(5), x2 = 2/sqrt(5)]];
190 block( [%rnum:0],
191 algsys([g526+((sqrt(a^2+4)-a)/2+a)*g525,
192 ((sqrt(a^2+4)-a)*g526)/2+g525],
193 [g525,g526]));
194 [[g525 = %r1,g526 = -((%r1*sqrt(a^2+4)+%r1*a)/2)]];
196 block( [%rnum:0],
197 algsys([g625+(a-(sqrt(a^2+4)+a)/2)*g624,
198 g624-((sqrt(a^2+4)+a)*g625)/2],
199 [g624,g625]));
200 [[g624 = %r1,g625 = (%r1*sqrt(a^2+4)-%r1*a)/2]];
202 /* Working example from bug #3150 */
203 algsys([x^2+y^2=1, x^2+(y-3)^2=9],[x,y]);
204 [[x = -sqrt(35)/6,y = 1/6],
205 [x = sqrt(35)/6,y = 1/6]];
207 /* Working example from bug #2736 */
208 algsys([x^2+y^2=1,(x-2)^2+(y-3)^2=10],[x,y]);
209 [[x = 1,y = 0],[x = -5/13,y = 12/13]];
211 /* bug #1106 only find 4 roots (of 10) but see below */
212 block( [
213 p1:-x*y^3+y^2+x^4-9*x/8,
214 p2:y^4-x^3*y-9*y/8+x^2],
215 algsys([p1,p2],[x,y]));
216 [[x = 1/2,y = 1],
217 [x = 9/8,y = 9/8],
218 [x = 1,y = 1/2],
219 [x = 0,y = 0]];
221 /* Above can be solved with
222 load(grobner)$
223 eqs:poly_reduced_grobner([p1,p2],[x,y])$
224 algsys(eqs,[x,y]);
225 */
226 algsys([-200704*y^8+437832*y^5-263633*y^2+53010*x,
227 4096*y^10-10440*y^7+7073*y^4-729*y],[x,y]);
228 [[x = 0, y = 0],
229 [x = 1, y = 1/2],
230 [x = 9/8, y = 9/8],
231 [x = 1/2, y = 1],
232 [x = (sqrt(3)*%i-1)/4, y = -((sqrt(3)*%i+1)/2)],
233 [x = -((sqrt(3)*%i+1)/4), y = (sqrt(3)*%i-1)/2],
234 [x = (3^(5/2)*%i-9)/16, y = -((3^(5/2)*%i+9)/16)],
235 [x = -((3^(5/2)*%i+9)/16), y = (3^(5/2)*%i-9)/16],
236 [x = (sqrt(3)*%i-1)/2, y = -((sqrt(3)*%i+1)/4)],
237 [x = -((sqrt(3)*%i+1)/2), y = (sqrt(3)*%i-1)/4]];
239 /* bug #1266 no solution found */
240 algsys([x^2+x+1=0,x^2*y+1=0],[x,y]);
241 [[x = (sqrt(3)*%i-1)/2, y = -((sqrt(3)*%i-1)/2)],
242 [x = -((sqrt(3)*%i+1)/2), y = (sqrt(3)*%i+1)/2]];
244 /* bug #1302 no solution found. Also bug #1469 */
245 algsys([(x-x0)^2+(y-y0)^2=r^2,(x-x1)^2+(y-y1)^2=r^2],[x,y]);
246 []; /* wrong, but at least it returns */
248 /* bug #1369 algsys hangs */
249 /*
250 eq1 : a*x + c*y + d*y^2/2 = 0;
251 eq2 : b*x + e*x^2 + f*y - g*y^3 = h;
252 algsys([eq1,eq2],[x,y]);
253 */
255 /* bug #1370 Can solve [F=0,F-G=0] but not [F,G]
256 there are some work arounds in the report
257 */
258 F:2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1;
259 2*x^3*y^2+2*a*x^2*y^2+3*x^2*y+2*a*x*y+x+1;
260 G:x^3*y^2+2*a*x^2*y^2+a^2*x*y^2+x^2*y+2*a*x*y+a^2*y+1;
261 x^3*y^2+2*a*x^2*y^2+a^2*x*y^2+x^2*y+2*a*x*y+a^2*y+1;
262 /* FIXME: Testsuite issue. Can't match form of solution
263 algsys([F=0,F-G=0],[x,y]);
264 [[x = -(((-a^3)+sqrt(a^2-4*a)*(a^2-2*a)+4*a^2)/(2*sqrt(a^2-4*a))),
265 y = -(1/(a*sqrt(a^2-4*a)))],
266 [x = -((a^3+sqrt(a^2-4*a)*(a^2-2*a)-4*a^2)/(2*sqrt(a^2-4*a))),
267 y = 1/(a*sqrt(a^2-4*a))]];
268 */
270 /* hangs
271 algsys([F=0,G=0],[x,y]);
272 [];
273 */
275 /* Can solve with
276 eqs:poly_reduced_grobner([F,G],[x,y])
277 algsys:(eqs,[x,y]);
278 */
279 /* FIXME: Testsuite issue. Can't match form of solution
280 algsys([(a^4-4*a^3)*y^2-1,(a^4-4*a^3)*y+2*x+a^2-2*a],[x,y]);
281 [[x = -((a*sqrt(a^2-4*a)+a^2-2*a)/2), y = sqrt(a^2-4*a)/(a^3-4*a^2)],
282 [x = (a*sqrt(a^2-4*a)-a^2+2*a)/2, y = -(sqrt(a^2-4*a)/(a^3-4*a^2))]];
283 */
285 /* Bug #1647 */
286 algsys([x^2+y^2-d^2, (x-h)^2+(y-k)^2-s^2], [x,y]);
287 []; /* no solution found, but at least it returns */
289 /*
290 Here is a workaround
291 load(grobner);
292 [x^2+y^2-d^2, (x-h)^2+(y-k)^2-s^2];
293 poly_reduced_grobner(%, [x,y]);
294 algsys(%, [x,y]);
295 */
296 /* FIXME: Testsuite issue. Can't match form of solution
297 algsys([2*k*y+2*h*x+s^2-k^2-h^2-d^2,
298 (4*k^2+4*h^2)*y^2+(k*(4*s^2-4*d^2)-4*k^3-4*h^2*k)*y+s^4-2*d^2*s^2
299 +h^2*((-2*s^2)+2*k^2-2*d^2)+k^2*(2*d^2-2*s^2)+k^4+h^4+d^4],[x,y]);
300 [[x = -(((k*sqrt((-s^4)+(2*k^2+2*h^2+2*d^2)*s^2-k^4+(2*d^2-2*h^2)*k^2-h^4
301 +2*d^2*h^2-d^4)
302 +h*s^2-(h*k^2)-h^3-(d^2*h))
303 /(2*k^2+2*h^2))),
304 y = (h*sqrt((-s^2)+2*d*s+k^2+h^2-d^2)*sqrt(s^2+2*d*s-k^2-h^2+d^2)
305 -k*s^2+k^3+(h^2+d^2)*k)
306 /(2*k^2+2*h^2)],
307 [x = (k*sqrt((-s^4)+(2*k^2+2*h^2+2*d^2)*s^2-k^4+(2*d^2-2*h^2)*k^2-h^4
308 +2*d^2*h^2-d^4)
309 -h*s^2+h*k^2+h^3+d^2*h)
310 /(2*k^2+2*h^2),
311 y = -(((h*sqrt((-s^2)+2*d*s+k^2+h^2-d^2)*sqrt(s^2+2*d*s-k^2-h^2+d^2)
312 +k*s^2-k^3+((-h^2)-d^2)*k)
313 /(2*k^2+2*h^2)))]];
314 */
316 /* SF bug #3210 ALGSYS does not return (regression)
317 * also #3212 ALGSYS throws error
318 *
319 * Note: %nrum is the counter for the %r variables introduced
320 * by algsys. Set it to 0 for reproducible results
321 */
322 block([%rnum:0], algsys([
323 4*y*(y^2-28),
324 8*y*(y^2-28)*z,
325 -16*y*(y^2-23)*z,
326 8*y*(y^2-10)*z,
327 32*y*(z-y),
328 32*y*(z+y),
329 4*y*(y^2-28),
330 8*y*(y^2-28)*z,
331 -16*y*(y^2-23)*z,
332 8*y*(y^2-10)*z,
333 y*((4*y^2-112)*z^2-8*y^4+(185-4*x^2)*y^2),
334 -2*y*((4*y^2-92)*z^2-2*y^4+(21-4*x^2)*y^2),
335 y*((4*y^2-40)*z^2+(1-4*x^2)*y^2)],
336 [x,y,z]));
337 [[x=%r1,y=0,z=%r2]];
339 block([%rnum:0], algsys([
340 -8*a*z,
341 -3*(y-1)*y*z,-3*z*(z+y-1)*(z+y),
342 2*(y-1)*y*z*(6*y*z-3*z-4*b*y^2-2*y^2+4*b*y+2*y),
343 -8*a*z,
344 -3*(y-1)*y*z,
345 -4*z*((w^2-2*v*w+v^2+2*b-2*a)*z+8*a*y-4*a),
346 -2*(y-1)*y*z*((6*y-3)*z+((-4*b)-2)*y^2+(4*b+2)*y),
347 4*z
348 *((2*w^2-4*v*w+2*v^2+2*b-1)*z^2
349 +(((-2*w^2)+4*v*w-2*v^2-8*b+8*a-2)*y+(2*w+2*v-2)*x-4*v*w+4*b-4*a+2)
350 *z-12*a*y^2+12*a*y-2*a),
351 -z*(((4*w^2-8*v*w+4*v^2+14)*y^2+(((-8*w)-8*v+8)*x+16*v*w-18)*y+4*x^2
352 -8*x+3)
353 *z^2
354 +(((-32*b)-16)*y^3+(48*b+24)*y^2+((-16*b)-8)*y)*z+(8*b-8*a+4)*y^4
355 +((-16*b)+16*a-8)*y^3+(8*b-8*a+4)*y^2),
356 -z*((4*w^2-8*v*w+4*v^2-1)*z^3+(((-16*w^2)+32*v*w-16*v^2-32*b)*y
357 +(16*w+16*v-16)*x-32*v*w+16*b+8)
358 *z^2
359 +((4*w^2-8*v*w+4*v^2+48*b-48*a+20)*y^2
360 +(((-8*w)-8*v+8)*x+16*v*w-48*b+48*a-24)
361 *y+4*x^2-8*x+8*b-8*a+4)
362 *z+32*a*y^3-48*a*y^2+16*a*y),
363 -2*z
364 *(((4*w^2-8*v*w+4*v^2+2)*y+((-4*w)-4*v+4)*x+8*v*w-3)*z^3
365 +(((-4*w^2)+8*v*w-4*v^2-24*b-8)*y^2
366 +((8*w+8*v-8)*x-16*v*w+24*b+12)*y-4*x^2+8*x-4*b-2)
367 *z^2+((16*b-16*a+8)*y^3+((-24*b)+24*a-12)*y^2+(8*b-8*a+4)*y)*z
368 +4*a*y^4-8*a*y^3+4*a*y^2)],
369 [v,w,x,y,z]));
370 [[v=%r1,w=%r2,x=%r3,y=%r4,z=0]];
372 /* Test problems distilled from share/contrib/diffequations tests
374 The odelin method for linear second order ODEs calls algys.
375 Changes to algsys capability change the odelin results.
376 A number of tests below were obtained by tracing algsys
377 while solving ODES with odelin (or contrib_ode).
378 */
380 /* Solved by odelin for Kamke ODE 2.265. New solution 2016-09 */
381 algsys(
382 [y = -1,
383 z = 1,
384 (-w^2)+2*v*w-v^2+5,2*((w+v-1)*x+4*w^2-10*v*w+4*v^2-22),
385 (-x^2)-(12*w+12*v-12)*x+((-2*w)-2*v+2)*x+2*x-26*w^2+80*v*w-26*v^2+161,
386 (-4*x^2)+((-8*w)-8*v+8)*x+8*x-4*w^2+24*v*w-4*v^2+48,
387 2*(3*x^2+(13*w+13*v-13)*x+(6*w+6*v-6)*x-6*x+22*w^2-82*v*w+22*v^2-157),
388 -2*((-6*x^2)-(12*w+12*v-12)*x+((-4*w)-4*v+4)*x+12*x-10*w^2+52*v*w
389 -10*v^2+100),
390 (-13*x^2)-(24*w+24*v-24)*x+((-26*w)-26*v+26)*x+26*x-41*w^2+182*v*w
391 -41*v^2+344],
392 [v,w,x,y,z]);
393 [[v = (sqrt(5)+5)/2, w = -((sqrt(5)-5)/2), x = 3, y = -1, z = 1],
394 [v = -((sqrt(5)-5)/2), w = (sqrt(5)+5)/2, x = 3, y = -1, z = 1],
395 [v = (sqrt(5)-3)/2, w = -((sqrt(5)+3)/2), x = -1, y = -1, z = 1],
396 [v = -((sqrt(5)+3)/2), w = (sqrt(5)-3)/2, x = -1, y = -1, z = 1],
397 [v = (sqrt(5)+1)/2, w = -((sqrt(5)-1)/2), x = 3, y = -1, z = 1],
398 [v = -((sqrt(5)-1)/2), w = (sqrt(5)+1)/2, x = 3, y = -1, z = 1],
399 [v = (sqrt(5)+1)/2, w = -((sqrt(5)-1)/2), x = -1, y = -1, z = 1],
400 [v = -((sqrt(5)-1)/2), w = (sqrt(5)+1)/2, x = -1, y = -1, z = 1]];
402 /* Solved by odelin for Murphy ODE 1.182. New solution 2016-10 */
403 /* FIXME: Testsuite issue. Can't match form of solution
404 algsys(
405 [x = -2*sqrt(-a1*a3),
406 4*sqrt(-a1*a3)*((-4*a1*a3*(2*y-z))-4*a0*a3*sqrt(-a1*a3)),
407 8*(-a1*a3)^(3/2)*(z^2-2*z-a2^2+2*a2)],
408 [x,y,z]);
409 [[x = -2*sqrt(-a1*a3),
410 y = -((a0*sqrt(-a1*a3)+a1*a2-2*a1)/(2*a1)),
411 z = 2-a2],
412 [x = -2*sqrt(-a1*a3),
413 y = -((a0*sqrt(-a1*a3)-a1*a2)/(2*a1)),
414 z = a2]];
415 */
417 /* Solved by odelin for Kamke ODE 1.105. New solution 2016-10 */
418 /* FIXME: Testsuite issue. Can't match form of solution
419 algsys(
420 [4*sqrt(-a*c)*((-4*a*c*(2*y-z))-4*a*sqrt(-a*c)*d),
421 8*(-a*c)^(3/2)*(z^2-2*z-b^2+2*b)],
422 [y,z]);
423 [[y = -((sqrt(-a*c)*d+(b-2)*c)/(2*c)), z = 2-b],
424 [y = -((sqrt(-a*c)*d-b*c)/(2*c)), z = b]];
425 */
427 /* Solved by odelin for Kamke ODE 1.291. New solution 2016-10 */
428 algsys([
429 -2*(18*y^2-36*x*y+18*x^2-5),
430 (72*y+72*x-72)*z+72*y^2-288*x*y+72*x^2+69,
431 (-36*z^2)-(144*y+144*x-144)*z+72*z-36*y^2+360*x*y-36*x^2-203,
432 72*z^2-((-72*y)-72*x+72)*z-144*z-144*x*y+159,(-36*z^2)+72*z-35],
433 [x,y,z]);
434 [[x = (sqrt(10)-3)/12, y = -((sqrt(10)+3)/12), z = 5/6],
435 [x = -((sqrt(10)+3)/12), y = (sqrt(10)-3)/12, z = 5/6],
436 [x = (sqrt(10)+13)/12, y = -((sqrt(10)-13)/12), z = 5/6],
437 [x = -((sqrt(10)-13)/12), y = (sqrt(10)+13)/12, z = 5/6],
438 [x = (sqrt(10)+15)/12, y = -((sqrt(10)-15)/12), z = 7/6],
439 [x = -((sqrt(10)-15)/12), y = (sqrt(10)+15)/12, z = 7/6],
440 [x = (sqrt(10)-1)/12, y = -((sqrt(10)+1)/12), z = 7/6],
441 [x = -((sqrt(10)+1)/12), y = (sqrt(10)-1)/12, z = 7/6]];
443 /* SF bug #2038
444 solve(eqs,[x,y]) gives correct solution
445 solve(subst(sqrt(3),q,eqs),[x,y]) yields []!
446 */
447 algsys(eqs:[y^2+x^2-(400-400/(2/q+1))^2=0,
448 (-y-400/(2/q+1)+400)^2+x^2-640000/(2/q+1)^2=0],[x,y]);
449 [[x = -(400*q*sqrt(-((q-2)/(q+2)))), y = -((400*q^2-800)/(q+2))],
450 [x = 400*q*sqrt(-((q-2)/(q+2))), y = -((400*q^2-800)/(q+2))]];
452 algsys(subst(q=sqrt(3),eqs),[x,y]);
453 [[x = 1200-800*sqrt(3),y = 400*sqrt(3)-800],
454 [x = 800*sqrt(3)-1200,y = 400*sqrt(3)-800]];
456 algsys(subst(q=1-sqrt(3),eqs),[x,y]);
457 [[x = -((25*2^(9/2))/3^(1/4)), y = 800/sqrt(3)],
458 [x = (25*2^(9/2))/3^(1/4), y = 800/sqrt(3)]];
460 /* SF bug #2736 (in comments)
461 Solving with p=1/(1-cos(1)) gave incorrect solution */
462 eqs:[x*y=p,x-y=3]$
463 [x*y = p,x-y = 3];
465 algsys(subst(p=1/(1-cos(1)),eqs),[x,y])$
466 [[x = -((sqrt(9*cos(1)^2-22*cos(1)+13)-3*cos(1)+3)/(2*cos(1)-2)),
467 y = -((sqrt(9*cos(1)^2-22*cos(1)+13)+3*cos(1)-3)/(2*cos(1)-2))],
468 [x = (sqrt(9*cos(1)^2-22*cos(1)+13)+3*cos(1)-3)/(2*cos(1)-2),
469 y = (sqrt(9*cos(1)^2-22*cos(1)+13)-3*cos(1)+3)/(2*cos(1)-2)]];
471 float(%)$
472 [[x = 3.603649840080336, y = 0.603649840080336],
473 [x = -0.603649840080336, y = -3.603649840080336]];
475 float(algsys(subst(p=float(1/(1-cos(1))),eqs),[x,y]))$
476 [[x = -0.6036497963666444,y = -3.603649796366644],
477 [x = 3.603649796366644, y = 0.6036497963666444]];
479 /* SF bug 2059 - variable order dependent
480 No solution found for certain solution variable orders
481 */
482 (eqs:[z-a=%i*b,z+a=c,a*z+3*(a-z)=13-12*%i,3*(z-a)+a*z=12*%i+13],
483 algsys(eqs,[a,b,c,z]));
484 [[a = (-2*%i)-3, b = 4, c = -6, z = 2*%i-3],
485 [a = 3-2*%i, b = 4, c = 6, z = 2*%i+3]];
487 /* Count the number of solutions for each permutation of variable */
488 map(length,
489 map(lambda([u],algsys(eqs,u)), listify(permutations([a,b,c,z]))));
490 [2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2];
492 /* SF bug #453 algsys fails in simple case */
493 algsys([v^2-2*v+u^2+10*u+1,v^2-10*v+u^2+2*u+1],[v,u]);
494 [[v = -((sqrt(34)-6)/2), u = (sqrt(2)*sqrt(17)-6)/2],
495 [v = (sqrt(34)+6)/2, u = -((sqrt(2)*sqrt(17)+6)/2)]];
497 algsys([2*v^2-12*v+1,v^2-2*v+u^2+10*u+1],[v,u]);
499 [[v = -((sqrt(2)*sqrt(17)-6)/2), u = -((sqrt(34)+14)/2)],
500 [v = (sqrt(2)*sqrt(17)+6)/2, u = (sqrt(34)-14)/2],
501 [v = (sqrt(2)*sqrt(17)+6)/2, u = -((sqrt(34)+6)/2)],
502 [v = -((sqrt(2)*sqrt(17)-6)/2), u = (sqrt(34)-6)/2]];
504 /* Solved by odelin for Kamke ODE 2.184. */
505 algsys(
506 [w^2-a^2*w,2*w*x-2*a^2*x,
507 w*(2*y+2*a*c+b^2)-3*a^2*y+x^2-2*a*b*x+A*a^2,
508 w*(2*z+4*b*c)-4*a^2*z+2*x*y-4*a*b*y+2*A*a*b,
509 x*(2*z+2*b*c)-6*a*b*z+y^2+((-2*a*c)-b^2)*y+3*c^2*w+2*A*a*c+A*b^2,
510 2*y*z+((-4*a*c)-2*b^2)*z+2*c^2*x+2*A*b*c,z^2-2*b*c*z+c^2*y+A*c^2],
511 [w,x,y,z]);
512 [[w = a^2,
513 x = -((a*sqrt((-4*a*c)+b^2-4*A)-3*a*b)/2),
514 y = -((b*sqrt((-4*a*c)+b^2-4*A)-2*a*c-b^2)/2),
515 z = -((c*sqrt((-4*a*c)+b^2-4*A)-b*c)/2)],
516 [w = a^2,
517 x = (a*sqrt((-4*a*c)+b^2-4*A)+3*a*b)/2,
518 y = (b*sqrt((-4*a*c)+b^2-4*A)+2*a*c+b^2)/2,
519 z = (c*sqrt((-4*a*c)+b^2-4*A)+b*c)/2]];
521 /* Solved by odelin for Kamke ODE 2.119. */
522 algsys([x^2-x-b^2-b,2*x*y-2*y,y^2],[x,y]);
523 [[x = -b, y = 0], [x = b+1, y = 0]];
525 /* Solved by odelin for Kamke ODE 2.331. */
526 algsys([x^2-16*x+48,2*x*y-32*y-192,y^2+64*y+64*x+192],[x,y]);
527 [[x = 12,y = -24],[x = 4,y = -8]];
529 /* Solved by odelin for Kamke ODE 2.382. */
530 algsys([
531 w^2-w,
532 2*w*x-2*x,w*(2*y+b^2+4*a*b+a^2)-3*y+x^2+(2*b+2*a)*x-c,
533 w*(2*z-4*a*b^2-4*a^2*b)-4*z+2*x*y+(4*b+4*a)*y+(2*b+2*a)*c,
534 x*(2*z-2*a*b^2-2*a^2*b)+(6*b+6*a)*z+y^2+((-b^2)-4*a*b-a^2)*y
535 +3*a^2*b^2*w+((-b^2)-4*a*b-a^2)*c,
536 2*y*z+((-2*b^2)-8*a*b-2*a^2)*z+2*a^2*b^2*x+(2*a*b^2+2*a^2*b)*c,
537 z^2+(2*a*b^2+2*a^2*b)*z+a^2*b^2*y-a^2*b^2*c],
538 [w,x,y,z]);
539 [[w = 1,
540 x = -((sqrt(4*c+b^2-(2*a*b)+a^2)+3*b+3*a)/2),
541 y = ((b+a)*sqrt(4*c+b^2-(2*a)*b+a^2)+b^2+4*a*b+a^2)/2,
542 z = -((a*b*sqrt(4*c+b^2-(2*a*b)+a^2)+a*b^2+a^2*b)/2)],
543 [w = 1,
544 x = (sqrt(4*c+b^2-(2*a*b)+a^2)-3*b-3*a)/2,
545 y = -(((b+a)*sqrt(4*c+b^2-2*a*b+a^2)-b^2-4*a*b-a^2)/2),
546 z = (a*b*sqrt(4*c+b^2-2*a*b+a^2)-a*b^2-a^2*b)/2]];
548 /* Solved by odelin for Kamke ODE 2.383. */
549 algsys([
550 w^2-4*w,
551 2*w*x-8*x,
552 w*(2*y+4*b^2+16*a*b+4*a^2)-12*y+x^2+(8*b+8*a)*x+b*(8*a*b^2-16*a^2*b+8*a^3)
553 +(4-4*a^2)*b^2+b^2*((-4*b^2)+8*a*b-4*a^2)
554 +(8*a^3-8*a)*b-4*a^4+4*a^2,
555 w*(2*z-16*a*b^2-16*a^2*b)-16*z+2*x*y+(16*b+16*a)*y
556 +b*((-16*a*b^3)+16*a^2*b^2+16*a^3*b-16*a^4)
557 +b^2*(8*b^3-8*a*b^2-8*a^2*b+8*a^3)+(8*a^2-8)*b^3
558 +(8*a-8*a^3)*b^2+(8*a^2-8*a^4)*b+8*a^5-8*a^3,
559 x*(2*z-8*a*b^2-8*a^2*b)+(24*b+24*a)*z+y^2+((-4*b^2)-16*a*b-4*a^2)*y
560 +12*a^2*b^2*w
561 +b*(8*a*b^4+16*a^2*b^3-48*a^3*b^2+16*a^4*b+8*a^5)
562 +(4-4*a^2)*b^4
563 +b^2*((-4*b^4)-8*a*b^3+24*a^2*b^2-8*a^3*b-4*a^4)
564 +(8*a-8*a^3)*b^3+(24*a^4-24*a^2)*b^2+(8*a^3-8*a^5)*b
565 -4*a^6+4*a^4,
566 2*y*z+((-8*b^2)-32*a*b-8*a^2)*z+8*a^2*b^2*x
567 +b*((-16*a^2*b^4)+16*a^3*b^3+16*a^4*b^2-16*a^5*b)
568 +b^2*(8*a*b^4-8*a^2*b^3-8*a^3*b^2+8*a^4*b)+(8*a^3-8*a)*b^4
569 +(8*a^2-8*a^4)*b^3+(8*a^3-8*a^5)*b^2+(8*a^6-8*a^4)*b,
570 z^2+(8*a*b^2+8*a^2*b)*z+4*a^2*b^2*y+b*(8*a^3*b^4-16*a^4*b^3+8*a^5*b^2)
571 +b^2*((-4*a^2*b^4)+8*a^3*b^3-4*a^4*b^2)+(4*a^2-4*a^4)*b^4
572 +(8*a^5-8*a^3)*b^3+(4*a^4-4*a^6)*b^2],
573 [w,x,y,z]);
574 [[w = 4,
575 x = 2*b^2+((-4*a)-6)*b+2*a^2-6*a,
576 y = (-2*b^3)+(2*a+2)*b^2+(2*a^2+8*a)*b-2*a^3+2*a^2,
577 z = 2*a*b^3+((-4*a^2)-2*a)*b^2+(2*a^3-2*a^2)*b],
578 [w = 4,
579 x = (-2*b^2)+(4*a-6)*b-2*a^2-6*a,
580 y = 2*b^3+(2-2*a)*b^2+(8*a-2*a^2)*b+2*a^3+2*a^2,
581 z = (-2*a*b^3)+(4*a^2-2*a)*b^2+((-2*a^3)-2*a^2)*b]];
583 /* Eugene L. Allgower, Kurt Georg and Rick Miranda, The Method of
584 Resultants for Computing Real Solutions of Polynomial Systems,
585 SIAM Journal on Numerical Analysis, Vol 29 No 3, Jun 1992, pp 831-844
586 */
588 /* Example 1, p840, qe 12. Robot arm coordinates.
589 Real solutions have y = +/- 1.6237
590 */
591 (p1:x^2*y^2+(2/q2)*x^2*y+(2/q2)*(1-b)*x*y^2+x^2+y^2-(2/q2)*(1+b)*x-(2/q2)*y+1,
592 p2:(1-b-q1)*x^2*y^2+(-1-b-q1)*x^2-4*x*y+(-1+b-q1)*y^2+(1+b-q1),
593 algsys(subst([b=1/2,q1=4/5,q2=2/5],[p1,p2]),[x,y]));
594 [[x = %i,y = -%i],
595 [x = -%i,y = %i],
596 [x = %i,y = -3*%i],
597 [x = -%i,y = 3*%i],
598 [x = (sqrt(319)+8)/17, y = -(sqrt(29)/sqrt(11))],
599 [x = -((sqrt(319)-8)/17), y = sqrt(29)/sqrt(11)]];
600 float([%[5],%[6]]);
601 [[x = 1.521210064675985, y = -1.623688281771977],
602 [x = -0.58003359408775, y = 1.623688281771977]];
604 /* Example 2, p 841. Intersection of a sphere with two paraboloids,
605 Two real solutions (not found by maxima 5.38.0)
606 y = z = (sqrt(5)-1)/2, x = +/- sqrt(z-y^2)
607 */
608 algsys([x^2+y^2+z^2-1,z-x^2-y^2,y-x^2-z^2],[x,y,z]);
609 [[x = sqrt(2)*%i, y = (sqrt(5)-1)/2, z = -((sqrt(5)+1)/2)],
610 [x = sqrt(2)*%i, y = -((sqrt(5)+1)/2), z = (sqrt(5)-1)/2],
611 [x = -sqrt(2)*%i, y = (sqrt(5)-1)/2, z = -((sqrt(5)+1)/2)],
612 [x = -sqrt(2)*%i, y = -((sqrt(5)+1)/2), z = (sqrt(5)-1)/2],
613 [x = sqrt((-sqrt(5))-2), y = -((sqrt(5)+1)/2), z = -((sqrt(5)+1)/2)],
614 [x = -sqrt((-sqrt(5))-2), y = -((sqrt(5)+1)/2), z = -((sqrt(5)+1)/2)],
615 [x = sqrt(sqrt(5)-2), y = (sqrt(5)-1)/2, z = (sqrt(5)-1)/2],
616 [x = -sqrt(sqrt(5)-2), y = (sqrt(5)-1)/2, z = (sqrt(5)-1)/2]];
618 /* Two examples from the maxima mailing list 2016-10-15 */
620 algsys(eqs:[x+2*y-z=6,2*x+y*z-z^2=-1,3*x-y+2*z^2=3],[x,y,z]);
621 [[x = (27*sqrt(13)*%i-41)/14,
622 y = -((17*sqrt(13)*%i-87)/14),
623 z = -((sqrt(13)*%i-7)/2)],
624 [x = -((27*sqrt(13)*%i+41)/14),
625 y = (17*sqrt(13)*%i+87)/14,
626 z = (sqrt(13)*%i+7)/2],
627 [x = 1, y = 2, z = -1]];
629 block([%rnum:0],
630 algsys(eqs:[b*c+a^2-1=0,b*d+a*b=0,c*d+a*c=0,d^2+b*c-1=0],[a,b,c,d]));
631 [[a = -sqrt(1-%r1*%r2), b = %r1, c = %r2, d = sqrt(1-%r1*%r2)],
632 [a = sqrt(1-%r3*%r4), b = %r3, c = %r4, d = -sqrt(1-%r3*%r4)],
633 [a = 1, b = %r5, c = 0, d = -1],
634 [a = -1, b = %r6, c = 0, d = 1],
635 [a = -1, b = 0, c = 0, d = -1],
636 [a = 1, b = 0, c = 0, d = 1]];
638 /* Example from Stavros Macrakis, maxima mailing list, 2017-07-31
639 Not solved by 5.38.0, solved by 5.39.0 */
640 block([%rnum:0],
641 algsys(eqs:[b*c+a^2,b*d+a*b,c*d+a*c,d^2+b*c],[a,b,c,d]));
642 [[a = -sqrt(-%r1*%r2), b = %r1, c = %r2, d = sqrt(-%r1*%r2)],
643 [a = sqrt(-%r3*%r4), b = %r3, c = %r4, d = -sqrt(-%r3*%r4)],
644 [a = 0, b = %r5, c = 0, d = 0]];
647 /* SF bug 3321 - apparent regression from 5.38.1
649 In fact, correct solution depends on the signs of constants g and h.
650 Version 5.38.1 returns wrong solution when g>0 and h<0.
651 Would be better if algsys asked for signs, but ...
652 */
654 eqs:[y-x=0, g*x*y-h=0, z+(x+1)/y-x-1=0];
655 [y-x=0, g*x*y-h=0, z+(x+1)/y-x-1=0];
657 (assume(h>0),0);
658 0;
659 s:algsys(eqs,[x,y,z]);
660 [[x = sqrt(h)/sqrt(g),
661 y = sqrt(h)/sqrt(g),
662 z = (h-g)/(sqrt(g)*sqrt(h))],
663 [x = -(sqrt(h)/sqrt(g)),
664 y = -(sqrt(h)/sqrt(g)),
665 z = -((h-g)/(sqrt(g)*sqrt(h)))]];
666 ratsimp(subst(s[1],eqs));
667 [0=0, 0=0, 0=0];
668 ratsimp(subst(s[2],eqs));
669 [0=0, 0=0, 0=0];
670 (forget(h>0),0);
671 0;
673 (assume(h<0),0);
674 0;
675 s:algsys(eqs,[x,y,z]);
676 [[x = -((%i*sqrt(-h))/sqrt(g)),
677 y = -((%i*sqrt(-h))/sqrt(g)),
678 z = -((sqrt(-h)*(%i*h-%i*g))/(sqrt(g)*h))],
679 [x = (%i*sqrt(-h))/sqrt(g),
680 y = (%i*sqrt(-h))/sqrt(g),
681 z = (sqrt(-h)*(%i*h-%i*g))/(sqrt(g)*h)]];
682 ratsimp(subst(s[1],eqs));
683 [0=0, 0=0, 0=0];
684 ratsimp(subst(s[2],eqs));
685 [0=0, 0=0, 0=0];
686 (forget(h<0),0);
687 0;
689 /* SF bug 3375. Failure to calculate eigenvalues
690 Solution is ugly, so just check it contains %r1 */
691 block([%rnum:0], s:algsys(eqs:[
692 y/10-(((-(sqrt(3)*%i)/2)-1/2)*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)
693 +((sqrt(3)*%i)/2-1/2)/(3*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)))*z,
694 x/20-(((-(sqrt(3)*%i)/2)-1/2)*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)
695 +((sqrt(3)*%i)/2-1/2)/(3*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)))*y,
696 50*z+20*y+((-((-(sqrt(3)*%i)/2)-1/2)*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3))
697 -((sqrt(3)*%i)/2-1/2)/(3*((sqrt(37)*%i)/(8*3^(3/2))+1/8)^(1/3)))*x],
698 [x,y,z]),0);
699 0;
700 /* Solution of homogeneous eqns must have free parameter */
701 freeof(%r1,s);
702 false;
704 /* SF bug 380. Missed solutions for algsys([a*b-c*d],[a,b,c,d]) */
705 block([%rnum:0],algsys([a*b-c*d],[a,b,c,d]));
706 [[a = (%r2*%r3)/%r1,b = %r1,c = %r2,d = %r3],
707 [a = %r4,b = 0,c = 0,d = %r5],
708 [a = %r6,b = 0,c = %r7,d = 0]];
710 block([%rnum:0],algsys([a*b*c-d*e*f],[a,b,c,d,e,f]));
711 [[a = (%r3*%r4*%r5)/(%r1*%r2),b = %r1,c = %r2,d = %r3,e = %r4,f = %r5],
712 [a = %r6,b = 0,c = %r7,d = 0,e = %r8,f = %r9],
713 [a = %r10,b = %r11,c = 0,d = 0,e = %r12,f = %r13],
714 [a = %r14,b = 0,c = %r15,d = %r16,e = 0,f = %r17],
715 [a = %r18,b = %r19,c = 0,d = %r20,e = 0,f = %r21],
716 [a = %r22,b = 0,c = %r23,d = %r24,e = %r25,f = 0],
717 [a = %r26,b = %r27,c = 0,d = %r28,e = %r29,f = 0]];
719 block([%rnum:0],algsys([a^2*b^2-c*d],[a,b,c,d]));
720 [[a = sqrt(%r2*%r3)/%r1,b = %r1,c = %r2,d = %r3],
721 [a = -sqrt(%r5*%r6)/%r4,b = %r4,c = %r5,d = %r6],
722 [a = %r7,b = 0,c = 0,d = %r8],
723 [a = %r9,b = 0,c = %r10,d = 0]];
725 block([%rnum:0],algsys([a*b^2*c-d*e],[a,b,c,d,e]));
726 [[a = (%r3*%r4)/(%r1^2*%r2),b = %r1,c = %r2,d = %r3,e = %r4],
727 [a = %r5,b = 0,c = %r6,d = 0,e = %r7],
728 [a = %r8,b = %r9,c = 0,d = 0,e = %r10],
729 [a = %r11,b = 0,c = %r12,d = %r13,e = 0],
730 [a = %r14,b = %r15,c = 0,d = %r16,e = 0]];
732 /* Reduced from rtest_to_poly_solve.mac
733 Third solution found after bug #380 fixed */
734 block([%rnum:0],
735 eqs:[a5*b1=0,a1*b2=r3,b2+a1=r2,a2*b2+a3*b1=r4,r4+a3*b2+a5=d5],
736 algsys(eqs,[a1,a2,a3,a5,b1,b2,r2,r3,r4]))$
737 [[a1=%r1,a2=%r2,a3=%r3,a5=d5+((-%r3)-%r2)*%r4,b1=0,b2=%r4,
738 r2=%r4+%r1,r3=%r1*%r4,r4=%r2*%r4],
739 [a1=%r5,a2=(d5-%r6*%r8-%r6*%r7)/%r8,a3=%r6,a5=0,b1=%r7,
740 b2=%r8,r2=%r8+%r5,r3=%r5*%r8,r4=d5-%r6*%r8],
741 [a1=%r9,a2=%r10,a3=d5/%r11,a5=0,b1=%r11,b2=0,r2=%r9,
742 r3=0,r4=d5]];
745 (kill(eqt,eqs,F,G,p,p1,p2,s),done)$
746 done;