Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

eliminate() #2720

Open
asmeurer opened this Issue Jan 1, 2014 · 11 comments

Comments

Projects
None yet
6 participants
@asmeurer
Copy link
Member

asmeurer commented Jan 1, 2014

We need some kind of eliminate function, like http://reference.wolfram.com/mathematica/ref/Eliminate.html. See also http://stackoverflow.com/q/20826969/161801 (and if someone, @smichr maybe, could answer that that would be great. I'm not sure what the best way is).

@anuraghota

This comment has been minimized.

Copy link

anuraghota commented Mar 6, 2014

Suppose for eliminate( x**2 +y**2 + z, x**2 - z , z) We can pass the 2 arguments as equation, solve the 2nd equation which gives something of form lets say { z : x**2} and replace z in first eqn with the value.

@anuraghota

This comment has been minimized.

Copy link

anuraghota commented Mar 6, 2014

e.g
1
2

And now replace it in first equation! Well there's some problem passing it as an equation .

@anuraghota

This comment has been minimized.

Copy link

anuraghota commented Mar 6, 2014

@asmeurer

This comment has been minimized.

Copy link
Member Author

asmeurer commented May 7, 2014

Another instance where I think this is what is wanted http://stackoverflow.com/q/23416592/161801

@jcrist jcrist added the enhancement label Apr 16, 2015

@Shekharrajak

This comment has been minimized.

Copy link
Member

Shekharrajak commented Jul 24, 2016

Any suggestions/idea, how it should be implemented ?

@asmeurer

This comment has been minimized.

Copy link
Member Author

asmeurer commented Jul 24, 2016

I think @smichr had some ideas. Do the answers on the stackoverflow question indicate an algorithm to do it?

@Shekharrajak

This comment has been minimized.

Copy link
Member

Shekharrajak commented Jul 25, 2016

Fine. I am trying it.

@sebant

This comment has been minimized.

Copy link

sebant commented Jul 29, 2016

I wrote one for my thermodynamics class a time ago. It doesn't follow any fancy strategy: it takes a list of equations and a list of symbols to eliminate, taking each symbol from the list finding a solution for it (from one of the equations) and updating the list by repeated substitution.
eliminate.zip

@Shekharrajak

This comment has been minimized.

Copy link
Member

Shekharrajak commented Jul 30, 2016

Thanks @sebant , but I don't want to use solve. I want to do this using subs, _invert .

@smichr

This comment has been minimized.

Copy link
Member

smichr commented Jul 1, 2017

I want to do this using subs, _invert .

This is my early attempt at this. I approached it from the point of view of trying to solve for some desired symbols rather than to eliminate some symbols:

def seek(eqs, do, sol=[], strict=True):
    from sympy.solvers.solvers import _invert as f
    from sympy.core.compatibility import ordered
    while do and eqs:
        for x in do:
            for e in eqs:
                i, d = f(e.lhs - e.rhs, x)
                if d != x:
                    continue
                break
            else:
                if strict:
                    assert None  # no eq could be solved for x
                continue
            sol.append((d, i))
            eqs.remove(e)
            break
        do.remove(x)
        if not strict:
            do.extend(i.free_symbols)
            do = list(ordered(do))
        for _ in range(len(eqs)):
            eqs[_] = eqs[_].xreplace({x: i})
    return sol

def focus(eqs, *syms, **kwargs):
    """Given Equality instances in ``eqs``, solve for symbols in
    ``syms`` and resolve as many of the free symbols in the solutions
    as possible. When ``evaluate=True`` a dictionary with keys being
    ``syms`` is returned, otherwise a list of identified symbols
    leading to the desired symbols is given.

    Examples
    ========
    >>> focus((Eq(a, b), Eq(b + 2, c)), a)
    {a: c - 2}
    >>> focus((Eq(a, b), Eq(b + 2, c)), a, evaluate=False)
    [(b, c - 2), (a, b)]
    """
    from sympy.solvers.solvers import _invert as f
    from sympy.core.compatibility import ordered
    evaluate = kwargs.get('evaluate', True)
    sol = []
    free = Tuple(*eqs).free_symbols
    do = set(syms) & free
    if not do:
        return sol
    eqs = list(eqs)
    seek(eqs, do, sol)
    assert not do
    for x, i in sol:
        do |= i.free_symbols
    do = list(ordered(do))  # make it canonical
    seek(eqs, do, sol, strict=False)
    if evaluate:
        while len(sol) > len(syms):
            x, s = sol.pop()
            for i in range(len(sol)):
                sol[i] = (sol[i][0], sol[i][1].xreplace({x: s}))
        for i in reversed(range(1, len(syms))):
            x, s = sol[i]
            for j in range(i):
                y, t = sol[j]
                sol[j] = y, f(y - t.xreplace({x: s}), y)[0]
    if evaluate:
        sol = dict(sol)
    else:
        sol = list(reversed(sol))
    return sol

Here are some results from a variety of places:

>>> x, y, z = var('x:z')
>>> eqs
[Eq(2*x + 3*y + 4*z - 1, 0), Eq(9*x + 8*y + 7*z - 2, 0)]
>>> focus(eqs, y, z, evaluate=1)
{z: x + 2/11, y: -2*x + 1/11}
>>> focus(eqs, y, z, evaluate=0)
[(y, -2*x + 1/11), (z, -x/2 - 3*y/4 + 1/4)]
>>> m1, m2 = symbols('m1, m2', real=True, positive=True)
>>> t0, t = symbols('t0 t', real=True, positive=True)
>>> p1, p2, v1, v2, a1, a2, f1, f2 = symbols('p1 p2 v1 v2 a1 a2 f1 f2', cls=Function)
>>> s = var('s')
>>> G = 6.67384e-11
>>>
>>> eq1 = Eq(f1(t), G * m1 * m2 * (p2(t) - p1(t)) / abs(p2(t) - p1(t)))
>>> eq2 = Eq(a1(t), a1(t0) + (1/m1) * Integral(f1(t), (t, t0, t)))
>>> eq3 = Eq(v1(t), v1(t0) + Integral(a1(t), (t, t0, t)))
>>> eq4 = Eq(p1(t), p1(t0) + Integral(v1(t), (t, t0, t)))
>>>
>>> eq5 = Eq(f2(t), G * m1 * m2 * (p1(t) - p2(t)) / abs(p2(t) - p1(t)))
>>> eq6 = Eq(a2(t), a2(t0) + (1/m2) * Integral(f2(t), (t, t0, t)))
>>> eq7 = Eq(v2(t), v2(t0) + Integral(a2(t), (t, t0, t)))
>>> eq8 = Eq(p2(t), p2(t0) + Integral(v2(t), (t, t0, t)))
>>> eqs = Tuple(eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8)
>>> eqs = eqs.subs((-p1(t) + p2(t))/Abs(p1(t) - p2(t)),s)
>>> eqs = eqs.subs((--p1(t) -+ p2(t))/Abs(p1(t) - p2(t)),-s)
>>> reps = zip(numbered_symbols(),eqs.atoms(Function))
>>> reps
[
(x0, v2(t)), 
(x1, p1(t0)), 
(x2, v2(t0)), 
(x3, p1(t)), 
(x4, p2(t)), 
(x5, v1(t)),
(x6, p2(t0)), 
(x7, f2(t)), 
(x8, a1(t)), 
(x9, v1(t0)), 
(x10, a2(t0)), 
(x11, a1(t0)), 
(x12, a2(t)), 
(x13, f1(t))]

>>> eqs = eqs.xreplace(dict([(v,k) for k, v in reps]))
>>> x3, x4 = [i for i,j in reps[3:5]]

>>> for k,v in focus(eqs,x3,x4,evaluate=0):(k,v)
...
(x5, x9 + Integral(x8, (t, t0, t)))
(x13, -1.0*x7)
(m1, -Integral(x13, (t, t0, t))/(x11 - x8))
(m2, 14983877347.9736*x13/(m1*s))
(x12, x10 + Integral(x7, (t, t0, t))/m2)
(x0, x2 + Integral(x12, (t, t0, t)))
(x4, x6 + Integral(x0, (t, t0, t)))
(x3, x1 + Integral(x5, (t, t0, t)))

>>> for k,v in focus(eqs,x3,x4).items():(k,v)
...
(x3, x1 + Integral(x9 + Integral(x8, (t, t0, t)), (t, t0, t)))
(x4, x6 + Integral(x2 + Integral(6.67384e-11*s*Integral(-1.0*x7, (t, t0, t))*Integral(x7, (t, t0, t))/(x7*(x11 - x8)) + x10, (t, t0, t)), (t, t0, t)))

# in terms of original functions

>>> for k, v in focus(eqs,x3,x4).items():(dict(reps)[k],v.xreplace(dict(reps)))
...
(p1(t), p1(t0) + Integral(v1(t0) + Integral(a1(t), (t, t0, t)), (t, t0, t)))
(p2(t), p2(t0) + Integral(v2(t0) + Integral(6.67384e-11*s*Integral(-1.0*f2(t), (t, t0, t))*Integral(f2(t), (t, t0, t))/((-a1(t) + a1(t0))*f2(t)) + a2(t0), (t, t0, t)), (t, t0, t)))
>>> e0, e1, e2, e3, e4, k, f, p, m = var('e:5 k f p m')
>>> for i in eqs: i
...
Eq(e1, k*(e0 - f + p) - m)
Eq(e2, k*(e1 - f + p) - m)
Eq(e3, k*(e2 - f + p) - m)
Eq(e4, -f + k*(e3 - f + p) - m)
>>> focus(eqs, e4)
{e4: -f + k*(-f + k*(-f + k*(-f + k*(e0 - f + p) - m + p) - m + p) - m + p) - m}

>>> for i in focus(eqs, e4, evaluate=False): i
...
(e1, k*(e0 - f + p) - m)
(e2, k*(e1 - f + p) - m)
(e3, k*(e2 - f + p) - m)
(e4, -f + k*(e3 - f + p) - m)
@asmeurer

This comment has been minimized.

Copy link
Member Author

asmeurer commented Dec 18, 2017

Is solving for subexpressions another instance where eliminate could be useful? See for instance https://stackoverflow.com/questions/47856322/sympy-solve-for-fraction. If you wanted to solve an expression for Y/X, could you set b = Y/X (where b is a dummy variable), and solve for b, eliminating Y and X? It seems like the systems solver should be able to be "tricked" into doing this, but I'm not quite sure how to best do it.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.