GSoC'20 - Week 3
In this week, work was done on two of the general case linear solvers type 2 and type 4. Namely:
- Type 2 solver: The PR which added the type 2 solver was merged this week. Most of the time taken on this PR was for
simplifying the solutions for some of the test cases but then it was concluded that instead of integrating
the part of the solution, we use the Integral function, thus an end user can use a
doit()function to integrate it if he/she chooses. This change led to successful passing of all the XFAIL test cases for this solver and thus this change was made. - Type 4 solver: The PR related to this solver was made this week and most of the base work related to adding a new solver like including the solver code, updating the matching function, adding test cases and determining the test cases that are slow was done. As was decided for type 2 solver, the Integral function is used instead of integrate function where an integration is done in the solution.
Linear, n equations, Order 1, Type 4: Non-constant coefficient non-homogeneous
The system is of the form below:
X' = A(t) * X + b(t)
Now, we will be using the same assumption that B, the antiderivative of A is commutative.
Hence:
We will use a solution similar to the solution for type 2 but instead of multiplying with exp(-A t), we will multiply with exp(-B(t)), i.e. exponential of the antiderivative.
exp(-B(t)) X' - A(t) exp(-B(t)) X = exp(-B(t)) f(t)
d/dt{ exp(-B(t)) X } = exp(-B(t)) f(t)
exp(-B(t)) X = int_x { exp(-B(t)) f(t) } + C
X = exp(B(t)) int_x { exp(-B(t)) f(t) } + exp(B(t)) C
The above solution is the one used by the solver to solve the systems of above type.
Lets look at an example to see the workings:
In [76]: eqs4
Out[76]:
⎡d d d ⎤
⎢──(f(x)) = x⋅(f(x) + g(x) + h(x)) + sin(x), ──(g(x)) = x⋅(f(x) + g(x) + h(x)) + sin(x), ──(h(x)) = x⋅(f(x) + g(x) + h(x)) + sin(x)⎥
⎣dx dx dx ⎦
Now, we have to obtain the matrices for solution:
In [64]: (A1, A0), b = linear_ode_to_matrix(eqs4, [f(x), g(x), h(x)], x, 1)
In [65]: A1
Out[65]:
⎡1 0 0⎤
⎢ ⎥
⎢0 1 0⎥
⎢ ⎥
⎣0 0 1⎦
In [66]: A0
Out[66]:
⎡x x x⎤
⎢ ⎥
⎢x x x⎥
⎢ ⎥
⎣x x x⎦
In [67]: b
Out[67]:
⎡sin(x)⎤
⎢ ⎥
⎢sin(x)⎥
⎢ ⎥
⎣sin(x)⎦
In [68]: B = integrate(A0, x)
In [69]: expB = B.exp()
Now, we have the commutative antiderivative B and the non-homogeneous term b. All we have to do now is get the solution using the formual described above:
In [70]: sol_vector = expB * (integrate(expB.inv() * b, x) + C)
In [75]: sol = [Eq(f, s) for f, s in zip([f(x), g(x), h(x)], sol_vector)]
The solution for this example is big and hence its not shown here. But, we can validate the solution:
In [77]: checksysodesol(eqs4, sol)
Out[77]: (True, [0, 0, 0])
These solutions were written by my mentor Oscar Benjamin, kudos to him. I have just implemented his methods for the project.
Future tasks
For the upcoming week, the target is to get the PR for type 4 solver merged and if possible, open a new PR for updating all linear solvers so that even an end user can directly access the solvers.