This week’s work involved adding solutions for more number of types of systems, two of which will be described in this week’s blog post.

Linear, n equations, Order 2, Type 2

The Order 2, Type 2 denotes system of the form given below:

X'' = f(t) * A * X
f(t) = 1/(a*t**2 + b*t + c)**2

Where, X is the vector of dependent variables, tis the dependent variable, Ais a constant matrix, and f(t) is the particular form of the function given below.

Lets understand the workflow with an example:

In [105]: eqs = [Eq(Derivative(f(t), (t, 2)), f(t)/t**4), Eq(Derivative(g(t), (t, 2)), d*g(t)/t**4)]                                                                                                        

In [106]: eqs                                                                                                                                                                                               
Out[106]: 
[Eq(Derivative(f(t), (t, 2)), f(t)/t**4),
 Eq(Derivative(g(t), (t, 2)), d*g(t)/t**4)]

In [107]: funcs = [f(t), g(t)]                                                                                                                                                                              

In [108]: (A2, A1, A0), b = linear_ode_to_matrix(eqs, [f(t), g(t)], t, 2) 

Our first check is to verify if A1 and b are 0 matrices or not.

In [111]: A1.is_zero_matrix, b.is_zero_matrix                                                                                                                                                               
Out[111]: (True, True)

The second check was to introduce a function that can detect such systems. For that, first step was to introduce a function named _factor_matrixwhich returns the factored matrix if the matrix is of the form: A(t) = f(t) * A_c where f(t) is a scalar function in t and A_c is a constant matrix else it returns None.

In [112]: _factor_matrix(A0, t)                                                                                                                                                                             
Out[112]: 
(t**(-4),
 Matrix([
 [1, 0],
 [0, d]]))

After we got the term f(t) and the factored matrix, the last and final check is to verify if the term f(t)is of the form 1/(a*t**2 + b*t+ c)**2. This is done internally by the function _is_second_order_type2. It uses the polynomial concepts to check if 1/f(t)is really a square of a quadratic function or not.

In [115]: _is_second_order_type2(A0, t)                                                                                                                                                                     
Out[115]: (True, t**2)

Now, after the check, this system is solved by using the substitutions: X(t) = g(t) * Y(tau), tau = integrate(f(t), t)where Y(tau) is a vector of new dependent variables with new independent variable tau and g(t) = 1/f(t)**(1/4)i.e. g(t) = sqrt(a*t**2 + b*t + c). This substitution transforms the original system into:

c1 = 1/4*b**2 - a*c
Y''(tau) = (A - I * c1) Y(tau)

where Iis the Identity matrix with same dimensions as A. Hence, now we have a constant coefficient matrix. This can be easily solved by introducing dummy variables and solving it separately. But after the solution is obtained, the solution will be in terms of tau. We have got the solution for Y(tau) and we need the solution of X(t). So, we use the old substitutions to get the answer:

X(t) = g(t) * Y(integrate(f(t), t)

Hence, we have our answer. The workings of these substitutions were internal and to keep the explanation simpler, no code was used to demonstrate it.

Euler Systems

Here, I have implemented the general solution for the cauhy euler systems, which are of the form:

t**n * X(n) * An + t**(n-1) *  X(n-1) * An-1 + ... + t**2 * X'' * A2 + t * X' * A1 + A0 * X = b(t)

To check if a higher order system is cauchy euler, one simply has to check that if all the coefficient matrices, divided by the appropriate term is constant or not. We do that by getting matrices from linear_ode_to_matrix.:

As, b = linear_ode_to_matrix(eqs, funcs, t, max_order)
is_euler = all(_matrix_is_constant(A/t**i) for i, A in enumerate(As))

This check is handled differently internally since the system input is in its canonical form, but the higher level logic of this check remains pretty much the same.

After we have done the check, its time to transform the above system using the substitutions X(t) = Y(tau), tau = ln(t) . This substitution transforms the system into constant coefficient system as described above:

X' = Y'/t
X'' = (Y'' - Y')/t**2
....

In this way, each non-constant term in the system becomes a constant one and the non-homogeneous term is transformed by doing the substitution b(tau) = b(t).subs(t, exp(tau)).

After solving for the constant coefficient system using the higher order to first order reduction, we need to obtain solution for X(t) from Y(tau) using:

X(t) = Y(ln(t))

This gives us the solution for the system.

Note

These cases are designed because the systems of above types, if reduced directly by using the dummy variables introduction method, don’t have a coefficient matrix which commutes with its own antiderivative.