GSoC'20 - Week 1 and Week 2
This blog contains the details of the work that I have done for week 1 and week 2. Since my college decided to take open book examinations online in the first week of June, I was busy in the first week with my exams. Hence, this post contains the work that I have done for both the weeks together.
Linear, n equations, Order 1, Type 2: Constant coefficient non-homogeneous solver
Lets look at the systems of this type:
X' = A*X + b(t)
This system looks very similar to the system described in type. Without going into too much detail, its easy to get the solution for the above system:
X' - AX = b(t)
exp(-A t) X' - exp(-A t) A X = exp(-A t)
Now, its important to note that A commutes with exp(-A t), so exp(-A t) A X = A exp(-A t) X:
exp(-A t) X' - A exp(-A t) X = exp(-A t) * b(t)
The lhs is just d/dt{ exp(-A t) X }. Hence,
d/dt{ exp(-A t) X } = exp(-A t) b(t)
Integrating on both the sides gives:
exp(-A t) X = integration( exp(-A t)*b(t), t ) + C
X = exp(A t)( integration( exp(-At)*b(t), t) + C )
Computing exponential of A t is expensive but we can use the jordan form method that we used for the solution for type 1.
exp(A t) = P * expJ * P.inv
exp(-A t) = P * expJ.inv * P.inv
Hence:
X = P * expJ * (integration( expJ.inv * P.inv * b(t), t ) + C)
Lets see the workings of this solution with the help of an example:
In [55]: eq1
Out[55]:
⎡d d ⎤
⎢──(f(x)) = f(x) + g(x) + 5, ──(g(x)) = -f(x) - g(x) + 7⎥
⎣dx dx ⎦
This is a constant coefficient non-homogeneous system. Lets get the coefficient matrices for finding the solution:
In [35]: (A1, A0), b = linear_ode_to_matrix(eq1, [f(x), g(x)], x, 1)
In [36]: A1
Out[36]:
⎡1 0⎤
⎢ ⎥
⎣0 1⎦
In [37]: A0
Out[37]:
⎡1 1 ⎤
⎢ ⎥
⎣-1 -1⎦
In [38]: b
Out[38]:
⎡5⎤
⎢ ⎥
⎣7⎦
We need only A0 and b for the solution. Now, we have to compute P and expJ.
Note: The method used here to compute the jordan normal form is not how its actually computed. There is a function created in the module before this GSoC project named matrix_exp_jordan_form specially designed to compute the exponential exp(A t) for faster computatiion given the fact that constant coeffcieint systems are very commonly in use.
In [40]: P, J = A0.jordan_form()
In [41]: P
Out[41]:
⎡1 1⎤
⎢ ⎥
⎣-1 0⎦
In [42]: J
Out[42]:
⎡0 1⎤
⎢ ⎥
⎣0 0⎦
In [43]: expJ = exp(J*x)
Now, we will use P, expJ and b to get the solution vector:
In [56]: sol_vector = P * expJ * ( integrate(expJ.inv() * P.inv() * b, x) + C)
In [57]: sol = [Eq(f, s) for f, s in zip([f(x), g(x)], sol_vector)]
In [59]: sol
Out[59]:
⎡ 2 2 ⎤
⎣f(x) = C₁ - 6⋅x - 7⋅x + (C₂ + 12⋅x)⋅(x + 1), g(x) = -C₁ + 6⋅x - x⋅(C₂ + 12⋅x) + 7⋅x
This is the solution for the system. Lets validate this solution:
In [58]: checksysodesol(eq1, sol)
Out[58]: (True, [0, 0])
Other tasks
Other than adding the above solver, some effort was made for simplifying the solutions by addition of two methods given below.
The simplification methods are:
-
simpsolmethod: Makes a polynomial like expression from the exponential terms which are factorized, simplifies the rational terms in all the coefficients of the made up polynomial. -
solsimpmethod: Simplifies the expression using powsimp, divides the expression into terms that are dependent and independent of the independent variable(t), carries out rational term simplification for the part of the expression independent oft, factorizes the power of the exponential terms of the expression dependent ontand adds both these simplified parts up.
These methods were checked throughly and analysis was also done for these methods. The link for the analysis is here.
In the upcoming week, the target is to complete the work for type 4 solver.