GSoC'20 - Week 6
This week, the work involved two PRs, one introduced dsolve_system that can solve system of ODEs. This function was designed in such a way so that to give a framework so that in future, when more solvers are added, it can be easily done so without changing the whole framework. Second PR involved the three tasks given by my mentor in past PRs, to clear the backlogs and make the framework better and more sensible.
PR 19695
This PR adds dsolve_system along with necessary helper functions.
Steps
The below steps shows how this PR handles the system of ODEs passed:
Step 1
Passes the system of ODEs through the matching function.
Step 2
If a match with one of the four types of system of ODEs introduced, then we directly use the
Step 3
_linear_ode_solver private function introduced, that takes in the match dictionary of a linear first order system of ODEs that has all the required details about the system and _linear_ode_solverinternally uses linodesolveto get the solution and return it.
If the system of ODEs is implicit or the match dictionary is None, then we pass the systems we got out of the canonical form(if the system is implicit, then we will have multiple systems from one single system, else we will just pass a list of single element containing the system which may or may not be solvable), we iterate throughout the systems in their canonical form and pass it to _component_solver. It is assumed that if a system is in its canonical form, then only one unique solution is possible: Picard–Lindelöf theorem
Step 4
Now, _component_solverfunction handles a system of ODEs, assuming that only one solution is possible. This function uses _component_division function and assumes that the division function divides the system of ODEs into strongly and weakly connected components, each component being a sub-system of the system of ODEs passed.
We will discuss more about how this work in the next week since for now, _component_division function is incomplete as it only returns a naive output which leads the _component_solver function and its helper function to assume that the system of ODEs cannot be divided into sub-systems and have to be solved as a whole. This is done so that we can work on the main thing of adding the framework as, if _component_division function was added in this PR itself, then this PR would have been extremely difficult to work with since it would become massive. Plus, now there would have been two things to track, the component division logic and the newly added framework. To avoid that, this function is made to output naive results for now. But, _component_solver function is made so that in future, when _component_division function is added, practically no changes have to be done in the framework or the _component_solver function.
The _component_solver function uses _weak_component_solver function to solve each weakly connected component and _weak_component_solveruses _strong_component_solver to solve each strongly connected component. Finally, all this is combined into one whole solution of the system passed after all the sub-systems are solved and this solution is returned back to dsolve_system function.
Step 5
After all the solutions have been obtained, it is time to prepare the solutions for final output. Now, linodesolvewas changed to return constants as dummies since otherwise, constants had to be tracked while solving different sub-systems, since if two sub-systems are solved separetly and they both have same constant symbols, then the solution becomes automatically incorrect. So, for each solution, these dummy symbols are replaced by constant symbols, if ics is provided(initial solution), then these constants are solved for and if doit is passed as True, then the unevaluated integrals are evaluated if any.
Types of systems of ODEs solved by dsolve_system
- Linear, First Order, Constant coefficient homogeneous system of ODEs
- Linear, First Order, Constant coefficient non-homogeneous system of ODEs
- Linear, First Order, non-constant coefficient homogeneous system of ODEs with coefficient matrix with commutative with its antiderivative
- Linear, First Order, non-constant coefficient non-homogeneous system of ODEs with coefficient matrix with commutative with its antiderivative
- Any implicit system which can be divided into system of ODEs which is of the above 4 forms
The internal working of the function has been demonstrated. Now, below are some examples of how an end-user can use the dsolve_system function to solve system of ODEs.
Examples
It is not mandatory to pass the dependent variables as a list or the symbol for independent variable as the function is made to automatically detect these variables. But a user can alterantively pass and get the same results.
In [34]: from sympy import symbols, Function, Eq
In [35]: from sympy.solvers.ode.systems import dsolve_system
In [36]: f, g = symbols("f g", cls=Function)
In [37]: system = [Eq(f(x).diff(x), f(x) + x), Eq(g(x).diff(x), f(x) + g(x))]
In [40]: dsolve_system(system)
Out[40]:
⎡⎡ ⌠ ⌠ ⌠ ⎤⎤
⎢⎢ x x ⎮ -x x x x ⎮ -x x ⎮ 2 -x ⎥⎥
⎢⎢f(x) = C₁⋅ℯ + ℯ ⋅⎮ x⋅ℯ dx, g(x) = C₁⋅x⋅ℯ + C₂⋅ℯ + x⋅ℯ ⋅⎮ x⋅ℯ dx + ℯ ⋅⎮ -x ⋅ℯ dx⎥⎥
⎣⎣ ⌡ ⌡ ⌡ ⎦⎦
In [41]: dsolve_system(system, funcs=[f(x), g(x)], t=x)
Out[41]:
⎡⎡ ⌠ ⌠ ⌠ ⎤⎤
⎢⎢ x x ⎮ -x x x x ⎮ -x x ⎮ 2 -x ⎥⎥
⎢⎢f(x) = C₁⋅ℯ + ℯ ⋅⎮ x⋅ℯ dx, g(x) = C₁⋅x⋅ℯ + C₂⋅ℯ + x⋅ℯ ⋅⎮ x⋅ℯ dx + ℯ ⋅⎮ -x ⋅ℯ dx⎥⎥
⎣⎣ ⌡ ⌡ ⌡ ⎦⎦
Now, as it is observable, this function passes a list of solutions, even when only one solution exists. If we pass an implicit system, which can be divided into multiple linear first order system of ODEs, then multiple solutions will be obtained.
In [42]: implicit_system = [Eq(f(x).diff(x)**2, g(x)**2), Eq(g(x).diff(x), f(x))]
In [43]: dsolve_system(implicit_system)
Out[43]:
⎡ ⎡ -x x -x x⎤⎤
⎣[f(x) = -C₁⋅sin(x) - C₂⋅cos(x), g(x) = C₁⋅cos(x) - C₂⋅sin(x)], ⎣f(x) = - C₁⋅ℯ + C₂⋅ℯ , g(x) = C₁⋅ℯ + C₂⋅ℯ ⎦⎦
As you can see, we have two solutions with us now. Both can be proved to satisfy the implicit system of ODEs. This has already been shown in the blog post about linodesolve hence the proof of its working is not being shown right now.
We can also pass doit=Trueas an option to evaluate the unevaluated integrals.
In [44]: system = [Eq(f(x).diff(x), f(x) + x), Eq(g(x).diff(x), f(x) + g(x))]
In [45]: dsolve_system(system, doit=True)
Out[45]:
⎡⎡ x x x 2 ⎤⎤
⎣⎣f(x) = C₁⋅ℯ - x - 1, g(x) = C₁⋅x⋅ℯ + C₂⋅ℯ + x + x⋅(-x - 1) + 2⋅x + 2⎦⎦
If, a user wants to solve the system of ODEs with initial conditions, then they can pass these intial values as a dictionary. For example, for the system used in these examples, if f(0) = 0 and g(0) = 1, then, we will pass this initial condition:
In [46]: dsolve_system(system, doit=True, ics={f(0): 0, g(0): 1})
Out[46]:
⎡⎡ x 2 x x ⎤⎤
⎣⎣f(x) = -x + ℯ - 1, g(x) = x + x⋅(-x - 1) + x⋅ℯ + 2⋅x - ℯ + 2⎦⎦
PR 19733
This PR extends constant_renumber and checkodesol and changes linear_ode_to_matrix to sync with canonical_odes function output. The PR is nearly done and it will be merged by today after adding some more test cases for coverage.
checkodesol
A checksysodesol function call was added in checkodesol which allows checkodesol to check the validity of the solutions of systems of ODEs. Below is an example demonstrating the same:
In [47]: system = [Eq(f(x).diff(x), f(x) + x), Eq(g(x).diff(x), f(x) + g(x))]
In [48]: sol = dsolve(system)
In [49]: checkodesol(system, sol)
Out[49]: (True, [0, 0])
constant_renumber
This function initially was not capable of handle system of dependent expressions and treated each expression separately. So, this would mean that if we pass a list of expressions:
In [3]: e1, e2, x, y = symbols("e1:3 x y")
In [4]: variables = [x, y]
In [7]: exprs = [e2*x, e1*x + e2*y]
And we want to replace e1 with C1 and e2with C2, then naturally, the solution should be either:
[C1*x, C2*x + C1*y] or [C2*x, C1*x + C2*y] but the answer is:
In [8]: constant_renumber(exprs, variables=variables)
Out[8]: [C1*x, C1*x + C2*y]
Which means that the function treats both the expressions separately. This is not desirable, hence this function was made to handle system of expressions.
So, with the updated function, the answer is:
In [5]: constant_renumber(exprs, variables=variables)
Out[5]: [C1*x, C1*y + C2*x]
linear_ode_to_matrix
This function initially returned coefficient matrices such that the system is of the form:
An * Xn + ... A1 * X1 + A0 * X = b
Where, Xn, ..., X1 are nth order, .., first order derivatives of the dependent variables vector, X is the vector of dependent variables and b is the non-homogeneous term. The function initially returned An, ..., A1, A0, bin the form given above.
Now, the function returns the same matrices but the system is now represented as:
An * Xn = An-1 * Xn-1 + ... A1 * X1 + A0 * X0 + b
Since, we want the coefficient matrices in the rhs when An is an identity matrix as our solvers assume these coefficient matrices as input to give appropriate solutions for the system of ODEs.
Below is an example demonstrating the correctness of the output of this function:
In [2]: f, g = symbols("f g", cls=Function)
In [3]: system = [Eq(f(x).diff(x), f(x) + x), Eq(g(x).diff(x), f(x) + g(x))]
In [4]: funcs = [f(x), g(x)]
In [5]: (A1, A0), b = linear_ode_to_matrix(system, funcs, x, 1)
In [6]: system_mat = Matrix([eq.lhs - eq.rhs for eq in system])
In [7]: system_mat
Out[7]:
Matrix([
[ -x - f(x) + Derivative(f(x), x)],
[-f(x) - g(x) + Derivative(g(x), x)]])
In [8]: X = Matrix(funcs)
In [9]: A1 * X.diff(x) - A0 * X - b == system_mat
Out[9]: True