Skip to main content

Week 11 and Week 12: Wrapping up the project

The GSoC is about to end. We have the suggested deadline tomorrow and the firm deadline next week. We are wrapping up the project now, cleaning up the code, writing the documentation and turning the project into a usable state. Since the new solvers cannot be called compelete I have to make sure that the next person working on it can pick it up easily. I plan come up with a post with the summary of the work that was done, that challenges we faced and thing that needs to done and how can they be tackled.

week 10: Radical equations

This week I worked on the solvers for the equations with radicals. Suppose you have to solve

$$ \sqrt{x} + x - 2 = 0 $$.

then you have to move the radical to the right hand side of the equation.

$$ x - 2 = - \sqrt{x} $$

and then square at both sides

$$ x^2 - 4x + 4 = x $$

Now the equation is a polynomial in \( x \) can be solved with usual polynomial solving methods. Note that squaring both sides produce some extra solutions and we will have to check all the solutions obtained against the original equation. If there are more than one radicals involved we may have to apply the method recursively. For example in solving \( \sqrt{2x + 9} - \sqrt{x + 1} - \sqrt{x + 4} = 0 \) the method will recurse twice.

To implement the method I tried a pattern matching approach. The squaring part is easy the tricky part is identifying which part to move to the right hand side. First I tried to match the expression with the form sqrt(p) + q but it failed even for case like 4*sqrt(x) + x - 2 because no pattern matched to it. I had to use a*sqrt(p) + q with the condition that the expression matched to a shouldn't be zero. Now I can simply move the expression matched with p and terms multiplicated with it to the RHS and square both the sides.

Notice that this method for solving sqrt equation can work with any radical equation, if it were cube root instead of sqrt I just had to cube both the sides. OK so how do I mathch that expression? I tried to pattern matching with assumptions on the wild symbols but it doesn't work. I tried to match with somthing like a*p**Rational(1, m) + q but this also didn't work out because Rational(1, m) raises TypeError no matter what the assumption on the variable are. There is a proposal for a new pattern matcher, I have not closely checked the details but it will be able to work with assumption. You can see the proposal on the wiki here and if it is implemented then things will be good but I can't wait for it. I had no other option to check term by term for rational power. Here's the implementation

def _has_rational_power(expr):
    a, p, q = Wild('a'), Wild('p'), Wild('q')
    pattern_match = expr.match(a*p**q)
    if pattern_match is None or pattern_match[a] is S.Zero:
        return (False, None)
    elif isinstance(pattern_match[q], Rational):
        if not pattern_match[q].q == S.One:
            return (True, pattern_match[q].q)

    if not isinstance(pattern_match[a], Pow) or isinstance(pattern_match[a], Mul):
        return (False, None)
    else:
        return _has_rational_power(pattern_match[a])

week 9

This week I moved back to college and my classes have restarted. This week I worked on a PR to allow infinitely indexed Range. See https://github.com/sympy/sympy/pull/7741. While doing this PR I discovered that you cannot monkey patch object to assign some attribute. And I want to mention that Sergey(one of my mentors) is damn good reviewer.

week 7

Week 7, 8

Fu Simplification

In the early part of the week 7 was thinking and working on design decisions. I want the code to be very modular so that it could be easily extended by people other than me. This came up in a meeting with Matthew and I want the solvers to be like the Fu simplification explained by Matthew in this Scipy Talk. The idea was that we can see solving an equation as series of transformations. If we have a lot of small transformations such that the input type is same as output type, and some notion of what makes a "better" output we can search though the list of transformations running one on top of other. I also posted about it on the mailing list which brought out some flaws in the preliminary design. The idea is pretty crude in the current stage and I'll have to look deeper into it, but not now.

I also discussed about the implementing a pattern dispached based solver suggested by F.B on the mailing list. But we decided that it will be better if we finish the equation solver by the current technique first.

Intersection with S.Reals

I decribed in the last post that one way to solve trigonometric equation is rewriting them in terms of \( exp \). But that is \( exp \) in the complex domain and the solution of \(exp(x) = a \) is \( \left\{i \left(2 \pi n + \arg{\left (a \right )}\right) + \log{\left (\left\lvert{a}\right\rvert \right )}\; |\; n \in \mathbb{Z}\right\} \). Hence we have to filter out real solutions from the obtained solutions. The filering is equivalent to the intersection of the solutions with the \( \mathbb{R} \) set. Suppose \( g(x) \) and \( h(x) \) are real valued functions and we have to perform $$ \mathbb{R} \cap \left\{g{\left (n \right )} + i h{\left (n \right )}\; |\; n \in \mathbb{Z}\right\} $$ then the answer will be simply $$ \left\{g{\left (n \right )}\; |\; n \in \left\{h{\left (n \right )} = 0\; |\; n \in \mathbb{Z}\right\}\right\} $$

Separate the real and imaginary parts and equate the imaginary to zero but the problem was with the assumptions on the symbols. For example while separating real and imaginary parts of the equation.

In[]: (n + I*n).as_imag_real()
Out[]:(re(n) - im(n), re(n) + im(n))

That is because n is by default complex, even in the Lambda(..).expr. I wrote some code to decide the assumption on the variable of imageset from the baseset. See PR 7694. There was another issue that needs to be resolved S.Integers.intersect(S.Reals) doesn't evaluate to S.Reals.

LambertW and Multivariate Solvers

The method to solve equation containing exp and log function is using the LambertW function. LambertW function is the inverse of \( x \exp(x) \). The function is multivariate function both for the real and complex domains and Sympy has only one branch implemented. This also leads us to loss of solutions. Aaron gave an example here. But I'm pretty unfamiliar with solving by LambertW and LambertW itself and it will take me some time to build an understanding of them. As an ad hoc solution I'm using the code in the _tsolve in the solvers.solvers module to do at least what the current solvers can do.

When the importing of _tsolve method was done. I started working on the multivariate solvers. Here's how the current multivariate solvers work:

Solving single multivariate equation

  1. count the number of free symbols in f - no of symbols for equation. If the equation has exactly one symbol which is not asked for then use solve_undetermined_coeffs, the solve_undetermined_coeffs find the values of the coefficient in a univariate polynomial such that it always equates to zero.

  2. Then for each symbol solve_linear is tried which tries to find a solution of that symbol in terms of constants or other symbols, the docstring says

    No simplification is done to f other than and mul=True expansion,
    so the solution will correspond strictly to a unique solution.
    

So we don't have to worry about loosing a solution. For every symbol it is checked if doesn't depend on previously solved symbols, if it does that solution is discarded.

  1. For the symbols for which the above method failed, the _solve function is called for the equation for that variable and as above if the solution contains a variable already solved then that solution is discarded.

System of equations in multiple variables

  • Try to convert the system of equations into a system of polynomial equation in variables

  • If all the equations are linear solve then using solve_linear_system, check the result and return it. If asked for particular solution solve using minsolve_linear_system

  • If the number of symbols is same as the size of the system solve the polynomial system using solve_poly_system. In case the system is over-determined All the free symbols intersection the variables asked for are calculated. Then for every subset of such symbols of length equal to that of the system, an attempt to solve the equations by solve_poly_system is made. Here if any of the solution depends on previously solved system the solution is discarded.

  • In the case there are failed equations:

    • For every know result:
    • Substitute every thing into the failed equation and see if the equation turns to zero. if it does accept the result otherwise put it in the bad_results group.
    • Then try to solve try to solve the failed equation using solve for each symbol.
    • If that solution depends on any other previously solved symbols discard it.
    • If it doesn't satisfy other equations, discard it.
    • Check if the solution doesn't set any denominator to zero, if it does discard that solution.
    • If it satisfies the above conditions substitute this value in know solutions and add it as a new result.

week 6

Solving Trigonometric Function (part II)

There is another technique to solve trigonometric function. Just as every trigonometric function can be written in term of \( \tan \) it can also be written in terms of \( \exp \).

$$ sin(x) = - \frac{i}{2} \left(e^{i x} - e^{- i x}\right) $$ $$ cos(x) = \frac{e^{i x}}{2} + \frac{1}{2} e^{- i x} $$ $$ tan(x) = \frac{i \left(- e^{i x} + e^{- i x}\right)}{e^{i x} + e^{- i x}} $$ $$ cot(x) = \frac{i \left(e^{i x} + e^{- i x}\right)}{e^{i x} - e^{- i x}} $$

So, solving a trigonometric equation is equivalent to solving a rational function in \( \exp \). Note: here the \( \exp \) is in complex domain and equation \( exp(x) = y \) has solution \( \left\{i \left(2 \pi n + \arg{\left (y \right )}\right) + \log{\left (\left\lvert{y}\right\rvert \right )}\; |\; n \in \mathbb{Z}\right\} \) when solved for \( x \).