GSoC'20 - Week 5
I was lucky enough to pass the first evaluation of GSoC’20. I would like to thank my mentors for their constant support and valuable feedback.
Aside from the first evaluations, this week, the main work done was combining the solvers introduced previously in this project and finally making a version of
dsolve_system, which will be able to take a system of ODEs as input and solve it appropriately using solvers added.
PR #19653
Work done in this PR is as follows:
- Combined the previously added solvers and combined it into a public function
linodesolve. - Added
linodesolve_typefunction, that can determine important information about the linear first-order system of ODEs, given the coefficient matrix, the independent variable, and the non-homogeneous term(if any). - Upgraded the
canonical_odesfunction and added support for the implicit system of ODEs.
This PR combined _linear_neq_order1_type1-4 solvers into linodesolve. A function that takes coefficient matrix, independent variable, and non-homogeneous term(if any), and matches the appropriate solution to solve for the type of system passed.
Let’s define a test system to test our linodesolve function on:
In [2]: f, g = symbols('f g', cls=Function)
...: eqs = [
Below are the things that are done in WEEK 4:
1. The PR related to Type 4 solver was merged. This PR was merged sooner than the other PRs since
this was the fourth time of adding a new solver.
2. It was observed that there are lots of repetitive lines in the four solvers added. Since it was
essential to work on them separately so that each of the solvers is tested out thoroughly and
worked upon individually. Now, these solvers are combined in the latest PR and this solver's
API is designed so that even an end-user can use the solver directly without going through the
matching function.
The target for WEEK 5 is to get the current PR merged as well as complete the evaluations properly.
...: Eq(f(x).diff(x), f(x) + y*g(x)),
...: Eq(g(x).diff(x), f(x)*y - g(x))]
...: funcs = [f(x), g(x)]
As you can see, the above system has the coefficient matrix defined below:
In [3]: M = Matrix([[1, y], [y, -1]])
Since we know that the system that we have is linear, order 1, constant-coefficient and non-homogeneous, the type of the system is “type1” and we can pass the appropriate arguments to get the solution for the system of ODEs:
In [4]: sol_vector = linodesolve(M, x, type="type1")
In [6]: sol_vector
Out[6]:
⎡ ________ ________ ⎤
⎢ ╱ 2 ╱ 2 ________ ________⎥
⎢ -x⋅╲╱ y + 1 x⋅╲╱ y + 1 ╱ 2 ╱ 2 ⎥
⎢ C₁⋅y⋅ℯ C₂⋅y⋅ℯ -x⋅╲╱ y + 1 x⋅╲╱ y + 1 ⎥
⎢- ──────────────────── + ───────────────────, C₁⋅ℯ + C₂⋅ℯ ⎥
⎢ ________ ________ ⎥
⎢ ╱ 2 ╱ 2 ⎥
⎣ ╲╱ y + 1 + 1 ╲╱ y + 1 - 1 ⎦
We can test the solution for the above system. The correctness of the solutions of the four types introduced is nothing to be worried since it is well tested in the file sympy/solvers/ode/tests/test_systems.py. But for the sake of testing this solution, the below is a method to test the correctness using checksysodesol
In [7]: sol = [Eq(func, s) for func, s in zip(funcs, sol_vector)]
In [8]: checksysodesol(eqs, sol)
Out[8]: (True, [0, 0])
(True, [0, 0]) means the solution for the system is correct.
linodesolve is designed to handle multiple cases. If the user doesn’t pass the type of the system that is to be solved, then with the help of the default value of the type keyword argument, the function will be able to determine the type of the system internally.
In [10]: sol_vector_without_type = linodesolve(M, x)
In [11]: sol_vector_with_type = linodesolve(M, x, type="type1")
In [13]: sol_vector_with_type == sol_vector_without_type
Out[13]: True
It is encouraged to pass the type of the system if the end-user knows the type of the system beforehand.
Since these types are defined internally, we created a function named linodesolve_type that takes the coefficient matrix, the independent variable, and the non-homogeneous term(if any) and returns the system information. In this information, the main thing is the type of the system, along with that it also returns the commutative antiderivative, if required to solve the system. of ODEs.
In [19]: f, g, h, k = symbols('f g h k', cls=Function)
...: eqs = [
...: Eq(f(x).diff(x), x*f(x) + x*y*g(x)),
...: Eq(g(x).diff(x), f(x)*y*x - x*g(x))]
...: funcs = [f(x), g(x)]
In [20]: M = Matrix([[x, x*y], [x*y, -1*x]])
In [21]: system_info = linodesolve_type(M, x)
In [22]: system_info
Out[22]:
{'type': 'type3',
'antiderivative': Matrix([
[ x**2/2, x**2*y/2],
[x**2*y/2, -x**2/2]])}
If the system is of Type 2 or Type 1, then the value for the key antiderivative is None.
In [25]: system_info = linodesolve_type(M/x, x)
In [26]: system_info
Out[26]: {'type': 'type1', 'antiderivative': None}
If a non-homogenous term is not passed(as seen above), or a zero vector is passed, then the system is assumed to be homogeneous. But, if a non-zero non-homogeneous term is passed, then the type of the system may be type2 or type4 depending on the fact that if the coefficient matrix is constant or not.
In [27]: linodesolve_type(M, x, b=None)['type']
Out[27]: 'type3'
In [28]: linodesolve_type(M, x, b=zeros(2, 1))['type']
Out[28]: 'type3'
In [29]: linodesolve_type(M, x, b=ones(2, 1))['type']
Out[29]: 'type4'
In [30]: linodesolve_type(M/x, x, b=ones(2, 1))['type']
Out[30]: 'type2'
Now, linodesolve cannot solve all the linear first-order system of ODEs. It can solve all constant coefficient systems but it can only solve a subset of the non-constant coefficient system of ODEs, namely, the systems whose coefficient matrix has a commutative antiderivative. If there is no commutative antiderivative for the coefficient matrix, then a NotImplementedError is raised by linodesolve_type(even by linodesolve if the system is non-constant, the antiderivative isn’t passed and the antiderivative isn’t commutative with its coefficient matrix).
In [31]: M = Matrix([[x, x*y], [y, -1]])
In [32]: linodesolve_type(M, x)
---------------------------------------------------------------------------
NotImplementedError Traceback (most recent call last)
...
NotImplementedError:
The system does not have a commutative antiderivative, it cannot be
solved by linodesolve.
There are system of implicit ODEs. Meaning, they can be solved for and divided into simpler systems and these systems may have a direct solution. We can get these solution systems from an implicit system using canonical_odes function.
The main objective of the canonical_odes when defined was to get a linear first order system of ODEs into the form:
` X’(t) = A(t) X(t) + b(t) `
where X(t) is a vector of dependent variables, t is the independent variable, A(t) is the coefficient matrix and b(t) is the non-homogeneous term.
This function was defined this way to easily get the coefficient matrix using linear_ode_to_matrix. We will demonstrate the use below:
In [33]: eqs = [Eq(f(x), f(x).diff(x)*2 + g(x)), Eq(0, -g(x) + g(x).diff(x)/12 + f(x))]
In [35]: canon_eqs = canonical_odes(eqs, funcs, x)
In [36]: canon_eqs
Out[36]:
⎡d f(x) g(x) d ⎤
⎢──(f(x)) = ──── - ────, ──(g(x)) = -12⋅f(x) + 12⋅g(x)⎥
⎣dx 2 2 dx ⎦
In [37]: eqs
Out[37]:
⎡ d ⎤
⎢ ──(g(x))⎥
⎢ d dx ⎥
⎢f(x) = g(x) + 2⋅──(f(x)), 0 = f(x) - g(x) + ────────⎥
⎣ dx 12 ⎦
In [38]: canon_eqs
Out[38]:
⎡d f(x) g(x) d ⎤
⎢──(f(x)) = ──── - ────, ──(g(x)) = -12⋅f(x) + 12⋅g(x)⎥
⎣dx 2 2 dx ⎦
In [39]: (A1, A0), b = linear_ode_to_matrix(eqs, funcs, x, 1)
In [40]: A1
Out[40]:
Matrix([
[-2, 0],
[ 0, -1/12]])
In [41]: A0
Out[41]:
Matrix([
[ 1, -1],
[-1, 1]])
In [42]: (A1, A0), b = linear_ode_to_matrix(canon_eqs, funcs, x, 1)
In [43]: A1
Out[43]:
Matrix([
[1, 0],
[0, 1]])
In [44]: A0
Out[44]:
Matrix([
[-1/2, 1/2],
[ 12, -12]])
We want the A1 matrix to be an identity matrix so that we can easily get the coefficient matrix in the desired form. The coefficient matrix to be passed if this system is to be solved is -A0. This is demonstrated below:
In [45]: (A1, A0), b = linear_ode_to_matrix(canon_eqs, funcs, x, 1)
In [46]: A = -A0
In [47]: sol_vector = linodesolve(A, x, b=b)
In [48]: sol = [Eq(func, s) for func, s in zip(funcs, sol_vector)]
In [49]: eqs
Out[49]:
⎡ d ⎤
⎢ ──(g(x))⎥
⎢ d dx ⎥
⎢f(x) = g(x) + 2⋅──(f(x)), 0 = f(x) - g(x) + ────────⎥
⎣ dx 12 ⎦
In [50]: checksysodesol(eqs, sol)
Out[50]: (True, [0, 0])
This is how dsolve used the linodesolve internally as well(with the help of a matching function).
As you can see, the main purpose of canonical_odes was to solve for the highest order derivative terms and give the solutions out. But, there can be multiple solutions as well. This functionality is added in canonical_odes.
In [51]: eqs = [Eq(f(x).diff(x)**2, g(x)**2), Eq(0, -g(x) + g(x).diff(x))]
In [52]: canon_eqs = canonical_odes(eqs, funcs, x)
In [53]: canon_eqs
Out[53]:
⎡⎡d d ⎤ ⎡d d ⎤⎤
⎢⎢──(f(x)) = -g(x), ──(g(x)) = g(x)⎥, ⎢──(f(x)) = g(x), ──(g(x)) = g(x)⎥⎥
⎣⎣dx dx ⎦ ⎣dx dx ⎦⎦
As you can see, we have two linear first order system of ODEs from one non-linear implicit system. This case is handled well by the dsolve function. It recognises the two distinct system of ODEs with the help of canonical_odes and solves both of those systems.
In [55]: sol = dsolve(eqs)
In [56]: for s in sol:
...: print(checksysodesol(eqs, s))
...:
(True, [0, 0])
(True, [0, 0])
It is easy enough to observe that both these solutions are correct. If the systems solved aren’t handled by linodesolve, then that system is passed to other solvers in dsolve.
This PR is merged and it was decided that two things needed to be rectified in the future:
- Adding
checksysodesolcall incheckodesolso that the user has to worry about only one solution checker. - Updating either
linear_ode_to_matrixorcanonical_odes, since, as of now, we have to take a negative sign ofA0to get the coefficient matrix that is to be used for the solution of the system of ODEs.
PR #19695
This PR adds the dsolve_system function that can solve for a system of ODEs passed. This PR is not merged yet so details about this PR will be discussed in the next week’s blog. The target is to get this PR merged this week and before moving on the other primary targets of the project, a PR addressing the targets assigned in the previous PR needs to be merged as well.