1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
|
/*
This routine is a modified version of linsolve written in maxima
Copyright (C) 1999 Dan Stanger
This library is free software; you can redistribute it and/or modify it
under the terms of the GNU Library General Public License as published
by the Free Software Foundation; either version 2 of the License, or (at
your option) any later version.
This library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Library General Public License for more details.
You should have received a copy of the GNU Library General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
Dan Stanger dan.stanger@ieee.org
please contact me for updates to this code
*/
/* modifies the global variable Solve_Inconsistent_Eqn_Nos if
Solve_Inconsistent_Error is false, LinsolveWarn and Linsolve_Params,
but acts as if Linsolve_Params is false */
/*
Added Solve_Inconsistent_Terms (December 2006, Andrej Vodopivec)
Solve_Inconsistent_Terms is not defined anywhere in solver package or
maxima. From how it is used I think it should be set in LinsolveM.
When triangularizing the coeff matrix sometimes rhs will be 0 and lhs
will be an expression in parameters. By setting Solve_Inconsistent_Terms
to be the list of these expressions, the user will be asked if they are
zero. If they are zero, the system can be solved (with the added
assumption that the expressions in Solve_Inconsistent_Terms are zero)
and if they are not zero, then the solver terminates - there are no
solutions.
With this change:
(%i38) Solver([x+y=p1-p2, x+3*y=p2-p3, 3*x-y=p1+p2+p3], [x,y], [p1,p2,p3]);
Is -p3+8*p2-4*p1 positive, negative, or zero? zero;
(%o38) [[x=-(-p3+4*p2-3*p1)/2,y=(-p3+2*p2-p1)/2]]
which I think is correct behaviour. Previously this triggered an error.
*/
define_variable( Solve_Inconsistent_Error, false, any)$
define_variable( Solve_Inconsistent_Eqn_Nos, false, any)$
define_variable( Solve_Inconsistent_Terms, [], list)$
LinsolveM(ees,vees):=block(
[ne:length(ees),nv:length(vees),nvees,ns,am,tam,r,
ic],
/* check to see if there are less equations than unknowns, if so,
solve for the first ones */
ns:min(ne,nv),
nvees:rest(vees,ns-nv),
/* construct the augmented coef matrix and triangularize it.
this assumes that all the non linear terms are on the rhs of the
equations (ees) */
am:augcoefmatrix(ees,nvees),
tam:triangularize(am),
/* search for inconsistent equations. this assumes that they are on
the bottom of the matrix. there may be a better way of doing
this but i am not sure of it */
/* copy the rhs, its the last column of the matrix, and whack it */
r:col(tam,ns+1),
tam:submatrix(tam,ns+1),
ic:[],
Solve_Inconsistent_Terms:[],
for i:ne thru 1 step -1 do block(
[az:true],
map(lambda([x],if x # 0 then az:false), tam[i]),
if az = true and r[i] # 0 then (
ic:endcons(i, ic),
Solve_Inconsistent_Terms : cons(ees[i], Solve_Inconsistent_Terms)
)
),
if ic # [] then ( /* equations are inconsistent */
if Solve_Inconsistent_Error = false then (
Solve_Inconsistent_Eqn_Nos:ic, return([inconsistant])
)
)
else (
if Solve_Inconsistent_Error = false then
Solve_Inconsistent_Eqn_Nos : []
else error("Inconsistent equations found!")
),
block(
[x:make_array('any,ns), l],
for i:ns thru 1 step -1 do
x[i-1]:(-r[i]-sum(tam[i,k]*x[k-1],k,i+1,ns))/tam[i,i],
l:map(first,listarray(x)),
if globalsolve = true then
/* this line does not work correctly */
map(lambda([x,y],x::y,globalsetq(x,y)),nvees,l)
else map(lambda([x,y],x=y),nvees,l)
)
)$
|
