undetermined_coeffs-like problem

[smichr]

I don't know how best to ask this question. If I have a equation with more unknowns that knowns and want to know the possible solutions and get something like this,

>>> solve(a*x+a+b*y-c-d, a, b)

[(0, (c + d)/y), (c/(x + 1), d/y), (d/(x + 1), c/y), ((c + d)/(x + 1), 0)]

Is there some significance of these 4 solutions? Are these solutions analogous to the case of `solve(a + b, a)` giving a = -b? The problem that I see is that if you write 

>>> solve(a*x+a+b*y-4, a, b)

[(0, 4/y), (4/(x + 1), 0)]

Those are only 2 of an infinite number of solutions since there are an infinite number of ways to partition up the 4. If we partition 4 as c + d (which then gives a solution that we obtained above) then substituting in different values of c and d gives a variety of solutions:

>>> ans=solve(a*x+a+b*y-c-d,a, b)

>>> [Tuple(*ai).subs({c:1,d:3}) for ai in ans]

[(0, 4/y), (1/(x + 1), 3/y), (3/(x + 1), 1/y), (4/(x + 1), 0)]

>>> [Tuple(*ai).subs({c:2,d:2}) for ai in ans]

[(0, 4/y), (2/(x + 1), 2/y), (2/(x + 1), 2/y), (4/(x + 1), 0)]

>>> [Tuple(*ai).subs({c:4,d:0}) for ai in ans]

[(0, 4/y), (4/(x + 1), 0), (0, 4/y), (4/(x + 1), 0)]

(Solving this type of equation arose in doing the multinomial_coefficients-without-expansion problem.)

/c

[Vinzent Steinberg]

Am Dienstag, 8. Mai 2012 10:44:04 UTC+2 schrieb smichr:

I don't know how best to ask this question. If I have a equation with more unknowns that knowns and want to know the possible solutions and get something like this,

>>> solve(a*x+a+b*y-c-d, a, b)

[(0, (c + d)/y), (c/(x + 1), d/y), (d/(x + 1), c/y), ((c + d)/(x + 1), 0)]

In this case, you have an underdetermined linear system, so you can solve this problem using linear algebra.

Suppose you want to solve

ax + a + by = z = c + d

for a and b. Written in matrix notation:

>>> A = Matrix([[x+1, y]])
>>> b = Matrix([z])

We want to find x, such that A*x == b.

Using Gaußian elimination, we obtain the reduced row echelon form of (A|b):

>>> A.row_join(b).rref()[0]
⎡     y      z  ⎤
⎢1  ─────  ─────⎥
⎣   x + 1  x + 1⎦

(Which is trivial in this case.)

Now we make the A part quadratic by inserting rows of the form (0...0, -1, 0...0) such that only 1 and -1 are on the diagonal:

>>> _.col_join(Matrix([[0, -1, 0]]))
⎡     y      z  ⎤
⎢1  ─────  ─────⎥
⎢   x + 1  x + 1⎥
⎢               ⎥
⎣0   -1      0  ⎦

This is all we need to obtain the solution space: The last column is a special solution to the inhomogeneous case (A*x == b). The linear combinations of columns which contain a -1 we added are all possible solutions to the homogeneous case (A*x == 0). So our solution is just the sum of the special inhomogeneous solution and the general homogeneous solution:

x = alpha * Matrix([y/(x + 1), -1]) + Matrix([z/(x +1), 0])

where alpha is any real number. (In this case our solution is a 1-dimensional vector space, in general it can have a higher dimension of course.)

Please note that you can obtain the homogeneous solution directly using A.nullspace(). I don't know whether there is already some function for the inhomogeneous solution.
 

Is there some significance of these 4 solutions? Are these solutions analogous to the case of `solve(a + b, a)` giving a = -b? The problem that I see is that if you write 

I get something else:

>>> solve(a*x + a + b*y - c - d, a, b)
⎡⎧   -b⋅y + c + d⎫⎤
⎢⎨a: ────────────⎬⎥
⎣⎩      x + 1    ⎭⎦

It does not really make sense to me to return exactly 4 solutions. However, you can use two special solutions to get the general one (constructing a line):

>>> x1 = Matrix([0, z/y])
>>> x2 = Matrix([z/(x + 1), 0])
>>> x2 - x1
⎡  z  ⎤
⎢─────⎥
⎢x + 1⎥
⎢     ⎥
⎢ -z  ⎥
⎢ ──  ⎥
⎣ y   ⎦
>>> _ / (z/y)
⎡  y  ⎤
⎢─────⎥
⎢x + 1⎥
⎢     ⎥
⎣ -1  ⎦

As you can see, alpha * (x2 - x1) + x2 is exactly the same result as above.
 

>>> solve(a*x+a+b*y-4, a, b)

[(0, 4/y), (4/(x + 1), 0)]

Those are only 2 of an infinite number of solutions since there are an infinite number of ways to partition up the 4. If we partition 4 as c + d (which then gives a solution that we obtained above) then substituting in different values of c and d gives a variety of solutions:

>>> ans=solve(a*x+a+b*y-c-d,a, b)

>>> [Tuple(*ai).subs({c:1,d:3}) for ai in ans]

[(0, 4/y), (1/(x + 1), 3/y), (3/(x + 1), 1/y), (4/(x + 1), 0)]

>>> [Tuple(*ai).subs({c:2,d:2}) for ai in ans]

[(0, 4/y), (2/(x + 1), 2/y), (2/(x + 1), 2/y), (4/(x + 1), 0)]

>>> [Tuple(*ai).subs({c:4,d:0}) for ai in ans]

[(0, 4/y), (4/(x + 1), 0), (0, 4/y), (4/(x + 1), 0)]

(Solving this type of equation arose in doing the multinomial_coefficients-without-expansion problem.)

If you always get linear systems, it is doable using the algorithm I described above.

Vinzent

[Chris Smith]

>
> I get something else:
>
>>>> solve(a*x + a + b*y - c - d, a, b)
> ⎡⎧   -b⋅y + c + d⎫⎤
> ⎢⎨a: ────────────⎬⎥
> ⎣⎩      x + 1    ⎭⎦
>
> It does not really make sense to me to return exactly 4 solutions. However,
> you can use two special solutions to get the general one (constructing a
> line):
>

Thanks, Vinzent. I've got to look at this in a unicode printout, but
let me at least respond to this part. The solve that I used is the
modified one in my multinomial branch. The 4 solutions I get are from
assigning all permutations of [c, d] to the a-part and b-part which
are being solved. e.g. the a-part is a*(x+1) and the b-part is b*y;
one solution is obtained by assigning -(c+d) to the a-part to give
a*(x+1)-(c+d)=0 and b*y=0 -> a = (c+d)/(x+1) and b=0, we can also
assign the -(c+d) to the b-part or the c to the a-part and d to the
b-part or vice versa. That's where the 4 solutions come from.

/c

[smichr]

Does anybody know if the following sorts of constraints can be solved to give the subsets of nonnegative integers that satisfy them:

i+j=2

k+l=4

m+n=3

i*k+l+m=3

j*k+l+n=8

The only thing I know to do is just check all possibilities:

def tallies(n, *syms):

    for p in partitions(n, len(syms)):

     take = sum([v for v in p.values()])

     vals = flatten([[k]*v for k, v in p.iteritems()])

     for s in subsets(syms, take):

      d = defaultdict(int)

      for i, si in enumerate(s):

          d[si] = vals[i]

      yield d

for dd in cartes(tallies(2,i,j), tallies(4,k,l), tallies(3,m,n)):

   d = defaultdict(int)

   for di in dd:

       d.update(di)

   if d[i]*d[k]+d[l]+d[m]==3 and d[j]*d[k]+d[l]+d[n]==8:

       print dict(d)

   

{k: 2, m: 1, j: 2, l: 2, n: 2, i: 0}

I suppose that it might be better to work backward from the more complicated constraints (like the i*k+l+m == 3).

Does anyone have other suggestions/algorithms?

/c

[Joachim Durchholz]

Am 10.05.2012 19:42, schrieb smichr:
> Does anybody know if the following sorts of constraints can be solved to
> give the subsets of nonnegative integers that satisfy them:
>
> i+j=2
> k+l=4
> m+n=3
> i*k+l+m=3
> j*k+l+n=8

This seems related to Diophantine equations.
Wikipedia has a bit to say about these - the article isn't deep but may
provide keywords for a websearch.

[Chris Smith]

Thanks for the reminder.

/c