Made Linearizer work with LagrangesMethod

This is the first step in getting the Linearizer Class to work with Lagranges Method. This required surprisingly few changes:

- Reformulation of general form to include lagrange multipliers.

- Change in `LagrangesMethod` interface to take in holonomic and nonholonomic constraints separately. Previously the holonomic constraints were time differentiated, combined with the nonholonomic constraints, and passed in via `coneqs`. 

Because this is a user facing change, this will have to handled with a deprecation warning. This has not been done yet, because I want to discuss the best method for handling this. Currently I've added `hol_coneqs` and `nonhol_coneqs` kwargs in the `__init__`. The original `coneqs` kwarg is still in there. If `coneqs` is present, it is used, otherwise it is assembled from the new kwargs. My thought is to add a deprecation warning when `coneqs` is used, and everything else should be good.

- Two tests exist using a minimal and nonminimal pendulum eom formulation. I plan to add more tests before merging.

**TODO:**

- [ ] Get the initial Linearizer class merged to sympy/master
- [ ] Documentation fixes for `LagrangesMethod` interface
- [ ] Documentation for `Linearizer` to update general form
- [ ] More tests. I think a rolling disc would be a good one.
- [ ] Figure out the issue of the Lagrange Multiplier operating point. See discussion [here](https://groups.google.com/forum/#!topic/pydy/j4vmsVizGvI) [1].
- [ ] Figure out the issue of expressing the dependent and independent variables [2]


**Reference [1]**
Basically, the linearization includes the values of the lagrange multipliers at the operating point. The question is how to go about doing this. The way I see it, there are 3 options:

1. Treat them like we treat all other operating point values, and have the user provide them in a dictionary. This is the current approach, and seems consistent.

2. Because the multipliers are not generally thought of as a system state, having them appear in the linearized form may be unexpected. They could be solved for inside `Linearizer`, and subbed out. This calculation could be computationally intensive however. Also, no other operating point conditions are solved for by `Linearizer`, so this may seem inconsistent.

3. Provide a `lam_value_at(operating_point)` method for the `LagrangesMethod` class, that solves for the multipliers at the provided operating point, and returns a dictionary of their values. The workflow would then be:

```
op_point = {dictionary of state: value}
lam_op = LM.lam_value_at(op_point)
op_point.update(lam_op)
linearizer = LM.to_linearizer()
M, A, B = linearizer.linearize(op_point)
```

**Reference [2]**
The issue here is that unlike `KanesMethod`, `LagrangesMethod` doesn't require knowledge of the dependent and independent variables in formulating the equations of motion. Currently I'm passing these into the `to_linearizer` method, but this is inconsistent with the other methods, and makes a general interface in the `linearize` function (not yet written) difficult. Perhaps `KanesMethod` could be rewritten to provide a formulation to `Linearizer` where the dependent and independent coordinates could be chosen later? This could be useful, and then both would have the same interface.