Skip to content

Commit

Permalink
Merge pull request #19341 from mijo2/const_coeff_nonhomogeneous
Browse files Browse the repository at this point in the history 
[GSOC] Constant coefficient non-homogeneous system of ODEs solver
  • Loading branch information
@oscarbenjamin
oscarbenjamin committed Jun 18, 2020
2 parents 8229061 + bf482b2 commit 25e8a79
Show file tree
Hide file tree
Showing 5 changed files with 450 additions and 369 deletions.
12 changes: 4 additions & 8 deletions doc/src/modules/solvers/ode.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -209,14 +209,6 @@ System of ODEs
These functions are intended for internal use by
:py:meth:`~sympy.solvers.ode.dsolve` for system of differential equations.

Linear, 2 equations, Order 1, Type 1
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. autofunction:: sympy.solvers.ode.ode::_linear_2eq_order1_type1

Linear, 2 equations, Order 1, Type 2
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. autofunction:: sympy.solvers.ode.ode::_linear_2eq_order1_type2

Linear, 2 equations, Order 1, Type 6
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. autofunction:: sympy.solvers.ode.ode::_linear_2eq_order1_type6
Expand Down Expand Up @@ -285,6 +277,10 @@ Linear, n equations, Order 1, Type 1
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. autofunction:: sympy.solvers.ode.systems::_linear_neq_order1_type1

Linear, n equations, Order 1, Type 2
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. autofunction:: sympy.solvers.ode.systems::_linear_neq_order1_type2

Linear, n equations, Order 1, Type 3
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. autofunction:: sympy.solvers.ode.systems::_linear_neq_order1_type3
Expand Down
213 changes: 9 additions & 204 deletions sympy/solvers/ode/ode.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -249,13 +249,13 @@
from sympy.core.symbol import Symbol, Wild, Dummy, symbols
from sympy.core.sympify import sympify

from sympy.logic.boolalg import (BooleanAtom, And, Not, BooleanTrue,
from sympy.logic.boolalg import (BooleanAtom, BooleanTrue,
BooleanFalse)
from sympy.functions import cos, cosh, exp, im, log, re, sin, sinh, sqrt, \
atan2, conjugate, Piecewise, cbrt, besselj, bessely, airyai, airybi
atan2, conjugate, cbrt, besselj, bessely, airyai, airybi
from sympy.functions.combinatorial.factorials import factorial
from sympy.integrals.integrals import Integral, integrate
from sympy.matrices import wronskian, Matrix
from sympy.matrices import wronskian
from sympy.polys import (Poly, RootOf, rootof, terms_gcd,
PolynomialError, lcm, roots, gcd)
from sympy.polys.polyroots import roots_quartic
Expand Down Expand Up @@ -2031,16 +2031,10 @@ def check_linear_2eq_order1(eq, func, func_coef):
# End of condition for type 6

if r['d1']!=0 or r['d2']!=0:
if not r['d1'].has(t) and not r['d2'].has(t):
if all(not r[k].has(t) for k in 'a1 a2 b1 b2 c1 c2'.split()):
# Equations for type 2 are Eq(a1*diff(x(t),t),b1*x(t)+c1*y(t)+d1) and Eq(a2*diff(y(t),t),b2*x(t)+c2*y(t)+d2)
return "type2"
else:
return None
return None
else:
if all(not r[k].has(t) for k in 'a1 a2 b1 b2 c1 c2'.split()):
# Equations for type 1 are Eq(a1*diff(x(t),t),b1*x(t)+c1*y(t)) and Eq(a2*diff(y(t),t),b2*x(t)+c2*y(t))
return "type1"
return None
else:
r['b1'] = r['b1']/r['a1'] ; r['b2'] = r['b2']/r['a2']
r['c1'] = r['c1']/r['a1'] ; r['c2'] = r['c2']/r['a2']
Expand All @@ -2050,6 +2044,7 @@ def check_linear_2eq_order1(eq, func, func_coef):
# Equations for type 7 are Eq(diff(x(t),t), f(t)*x(t) + g(t)*y(t)) and Eq(diff(y(t),t), h(t)*x(t) + p(t)*y(t))
return "type7"


def check_linear_2eq_order2(eq, func, func_coef):
x = func[0].func
y = func[1].func
Expand Down Expand Up @@ -6517,204 +6512,12 @@ def sysode_linear_2eq_order1(match_):
raise NotImplementedError("Only homogeneous problems are supported" +
" (and constant inhomogeneity)")

if match_['type_of_equation'] == 'type2':
gsol = _linear_2eq_order1_type1(x, y, t, r, eq)
psol = _linear_2eq_order1_type2(x, y, t, r, eq)
sol = [Eq(x(t), gsol[0].rhs+psol[0]), Eq(y(t), gsol[1].rhs+psol[1])]
if match_['type_of_equation'] == 'type6':
sol = _linear_2eq_order1_type6(x, y, t, r, eq)
if match_['type_of_equation'] == 'type7':
sol = _linear_2eq_order1_type7(x, y, t, r, eq)
return sol

# To remove Linear 2 Eq, Order 1, Type 1 when
# Linear 2 Eq, Order 1, Type 2 is removed.
def _linear_2eq_order1_type1(x, y, t, r, eq):
r"""
It is classified under system of two linear homogeneous first-order constant-coefficient
ordinary differential equations.
The equations which come under this type are
.. math:: x' = ax + by,
.. math:: y' = cx + dy
The characteristics equation is written as
.. math:: \lambda^{2} + (a+d) \lambda + ad - bc = 0
and its discriminant is `D = (a-d)^{2} + 4bc`. There are several cases
1. Case when `ad - bc \neq 0`. The origin of coordinates, `x = y = 0`,
is the only stationary point; it is
- a node if `D = 0`
- a node if `D > 0` and `ad - bc > 0`
- a saddle if `D > 0` and `ad - bc < 0`
- a focus if `D < 0` and `a + d \neq 0`
- a centre if `D < 0` and `a + d \neq 0`.
1.1. If `D > 0`. The characteristic equation has two distinct real roots
`\lambda_1` and `\lambda_ 2` . The general solution of the system in question is expressed as
.. math:: x = C_1 b e^{\lambda_1 t} + C_2 b e^{\lambda_2 t}
.. math:: y = C_1 (\lambda_1 - a) e^{\lambda_1 t} + C_2 (\lambda_2 - a) e^{\lambda_2 t}
where `C_1` and `C_2` being arbitrary constants
1.2. If `D < 0`. The characteristics equation has two conjugate
roots, `\lambda_1 = \sigma + i \beta` and `\lambda_2 = \sigma - i \beta`.
The general solution of the system is given by
.. math:: x = b e^{\sigma t} (C_1 \sin(\beta t) + C_2 \cos(\beta t))
.. math:: y = e^{\sigma t} ([(\sigma - a) C_1 - \beta C_2] \sin(\beta t) + [\beta C_1 + (\sigma - a) C_2 \cos(\beta t)])
1.3. If `D = 0` and `a \neq d`. The characteristic equation has
two equal roots, `\lambda_1 = \lambda_2`. The general solution of the system is written as
.. math:: x = 2b (C_1 + \frac{C_2}{a-d} + C_2 t) e^{\frac{a+d}{2} t}
.. math:: y = [(d - a) C_1 + C_2 + (d - a) C_2 t] e^{\frac{a+d}{2} t}
1.4. If `D = 0` and `a = d \neq 0` and `b = 0`
.. math:: x = C_1 e^{a t} , y = (c C_1 t + C_2) e^{a t}
1.5. If `D = 0` and `a = d \neq 0` and `c = 0`
.. math:: x = (b C_1 t + C_2) e^{a t} , y = C_1 e^{a t}
2. Case when `ad - bc = 0` and `a^{2} + b^{2} > 0`. The whole straight
line `ax + by = 0` consists of singular points. The original system of differential
equations can be rewritten as
.. math:: x' = ax + by , y' = k (ax + by)
2.1 If `a + bk \neq 0`, solution will be
.. math:: x = b C_1 + C_2 e^{(a + bk) t} , y = -a C_1 + k C_2 e^{(a + bk) t}
2.2 If `a + bk = 0`, solution will be
.. math:: x = C_1 (bk t - 1) + b C_2 t , y = k^{2} b C_1 t + (b k^{2} t + 1) C_2
"""

C1, C2 = get_numbered_constants(eq, num=2)
a, b, c, d = r['a'], r['b'], r['c'], r['d']
real_coeff = all(v.is_real for v in (a, b, c, d))
D = (a - d)**2 + 4*b*c
l1 = (a + d + sqrt(D))/2
l2 = (a + d - sqrt(D))/2
equal_roots = Eq(D, 0).expand()
gsol1, gsol2 = [], []

# Solutions have exponential form if either D > 0 with real coefficients
# or D != 0 with complex coefficients. Eigenvalues are distinct.
# For each eigenvalue lam, pick an eigenvector, making sure we don't get (0, 0)
# The candidates are (b, lam-a) and (lam-d, c).
exponential_form = D > 0 if real_coeff else Not(equal_roots)
bad_ab_vector1 = And(Eq(b, 0), Eq(l1, a))
bad_ab_vector2 = And(Eq(b, 0), Eq(l2, a))
vector1 = Matrix((Piecewise((l1 - d, bad_ab_vector1), (b, True)),
Piecewise((c, bad_ab_vector1), (l1 - a, True))))
vector2 = Matrix((Piecewise((l2 - d, bad_ab_vector2), (b, True)),
Piecewise((c, bad_ab_vector2), (l2 - a, True))))
sol_vector = C1*exp(l1*t)*vector1 + C2*exp(l2*t)*vector2
gsol1.append((sol_vector[0], exponential_form))
gsol2.append((sol_vector[1], exponential_form))

# Solutions have trigonometric form for real coefficients with D < 0
# Both b and c are nonzero in this case, so (b, lam-a) is an eigenvector
# It splits into real/imag parts as (b, sigma-a) and (0, beta). Then
# multiply it by C1(cos(beta*t) + I*C2*sin(beta*t)) and separate real/imag
trigonometric_form = D < 0 if real_coeff else False
sigma = re(l1)
if im(l1).is_positive:
beta = im(l1)
else:
beta = im(l2)
vector1 = Matrix((b, sigma - a))
vector2 = Matrix((0, beta))
sol_vector = exp(sigma*t) * (C1*(cos(beta*t)*vector1 - sin(beta*t)*vector2) + \
C2*(sin(beta*t)*vector1 + cos(beta*t)*vector2))
gsol1.append((sol_vector[0], trigonometric_form))
gsol2.append((sol_vector[1], trigonometric_form))

# Final case is D == 0, a single eigenvalue. If the eigenspace is 2-dimensional
# then we have a scalar matrix, deal with this case first.
scalar_matrix = And(Eq(a, d), Eq(b, 0), Eq(c, 0))

vector1 = Matrix((S.One, S.Zero))
vector2 = Matrix((S.Zero, S.One))
sol_vector = exp(l1*t) * (C1*vector1 + C2*vector2)
gsol1.append((sol_vector[0], scalar_matrix))
gsol2.append((sol_vector[1], scalar_matrix))

# Have one eigenvector. Get a generalized eigenvector from (A-lam)*vector2 = vector1
vector1 = Matrix((Piecewise((l1 - d, bad_ab_vector1), (b, True)),
Piecewise((c, bad_ab_vector1), (l1 - a, True))))
vector2 = Matrix((Piecewise((S.One, bad_ab_vector1), (S.Zero, Eq(a, l1)),
(b/(a - l1), True)),
Piecewise((S.Zero, bad_ab_vector1), (S.One, Eq(a, l1)),
(S.Zero, True))))
sol_vector = exp(l1*t) * (C1*vector1 + C2*(vector2 + t*vector1))
gsol1.append((sol_vector[0], equal_roots))
gsol2.append((sol_vector[1], equal_roots))
return [Eq(x(t), Piecewise(*gsol1)), Eq(y(t), Piecewise(*gsol2))]


def _linear_2eq_order1_type2(x, y, t, r, eq):
r"""
The equations of this type are
.. math:: x' = ax + by + k1 , y' = cx + dy + k2
The general solution of this system is given by sum of its particular solution and the
general solution of the corresponding homogeneous system is obtained from type1.
1. When `ad - bc \neq 0`. The particular solution will be
`x = x_0` and `y = y_0` where `x_0` and `y_0` are determined by solving linear system of equations
.. math:: a x_0 + b y_0 + k1 = 0 , c x_0 + d y_0 + k2 = 0
2. When `ad - bc = 0` and `a^{2} + b^{2} > 0`. In this case, the system of equation becomes
.. math:: x' = ax + by + k_1 , y' = k (ax + by) + k_2
2.1 If `\sigma = a + bk \neq 0`, particular solution is given by
.. math:: x = b \sigma^{-1} (c_1 k - c_2) t - \sigma^{-2} (a c_1 + b c_2)
.. math:: y = kx + (c_2 - c_1 k) t
2.2 If `\sigma = a + bk = 0`, particular solution is given by
.. math:: x = \frac{1}{2} b (c_2 - c_1 k) t^{2} + c_1 t
.. math:: y = kx + (c_2 - c_1 k) t
"""
r['k1'] = -r['k1']; r['k2'] = -r['k2']
if (r['a']*r['d'] - r['b']*r['c']) != 0:
x0, y0 = symbols('x0, y0', cls=Dummy)
sol = solve((r['a']*x0+r['b']*y0+r['k1'], r['c']*x0+r['d']*y0+r['k2']), x0, y0)
psol = [sol[x0], sol[y0]]
elif (r['a']*r['d'] - r['b']*r['c']) == 0 and (r['a']**2+r['b']**2) > 0:
k = r['c']/r['a']
sigma = r['a'] + r['b']*k
if sigma != 0:
sol1 = r['b']*sigma**-1*(r['k1']*k-r['k2'])*t - sigma**-2*(r['a']*r['k1']+r['b']*r['k2'])
sol2 = k*sol1 + (r['k2']-r['k1']*k)*t
else:
# FIXME: a previous typo fix shows this is not covered by tests
sol1 = r['b']*(r['k2']-r['k1']*k)*t**2 + r['k1']*t
sol2 = k*sol1 + (r['k2']-r['k1']*k)*t
psol = [sol1, sol2]
return psol

def _linear_2eq_order1_type6(x, y, t, r, eq):
r"""
The equations of this type of ode are .
Expand Down Expand Up @@ -7478,10 +7281,12 @@ def _linear_2eq_order2_type11(x, y, t, r, eq):

def sysode_linear_neq_order1(match):
from sympy.solvers.ode.systems import (_linear_neq_order1_type1,
_linear_neq_order1_type3)
_linear_neq_order1_type3, _linear_neq_order1_type2)

if match['type_of_equation'] == 'type1':
sol = _linear_neq_order1_type1(match)
elif match['type_of_equation'] == 'type2':
sol = _linear_neq_order1_type2(match)
elif match['type_of_equation'] == 'type3':
sol = _linear_neq_order1_type3(match)
return sol
Expand Down

0 comments on commit 25e8a79

Please sign in to comment.