Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[GSOC] Constant coefficient non-homogeneous system of ODEs solver #19341

Merged
merged 43 commits into from Jun 18, 2020

Conversation

mijo2
Copy link
Contributor

@mijo2 mijo2 commented May 17, 2020

References to other Issues or PRs

Fixes #9244
Fixes #8859
Fixes #8567
Fixes #19150

Brief description of what is fixed or changed

Added linear first-order constant-coefficient non-homogeneous system of ODEs solver.

Other comments

The SymPy's ode module will be able to solve a number of linear first-order constant coefficient non-homogenous systems of ODEs

For example:

>>> from sympy import *
>>> from sympy.abc import *
>>> f, g = symbols("f g", cls=Function)
>>> eqs = [Eq(f(t).diff(t), f(t) + 2*g(t) + t), Eq(g(t).diff(t), 7*g(t) + t**2)]
>>> dsolve(eqs)
[Eq(f(t), C1*exp(t) + C2*exp(7*t)/3 + 2*t**2/7 - 17*t/49 - 115/343), Eq(g(t), C2*exp(7*t) - t**2/7 - 2*t/49 - 2/343)]

Release Notes

  • solvers
    • In dsolve there is no a general solver that can solve systems of constant coefficient non-homogeneous first order ODEs of any size in terms of integrals.

Added the constant coefficient non-homogeneous solver.
@sympy-bot
Copy link

sympy-bot commented May 17, 2020

Hi, I am the SymPy bot (v160). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

  • solvers
    • In dsolve there is no a general solver that can solve systems of constant coefficient non-homogeneous first order ODEs of any size in terms of integrals. (#19341 by @mijo2)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.7.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it.

Click here to see the pull request description that was parsed.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
Fixes #9244
Fixes #8859
Fixes #8567
Fixes #19150


#### Brief description of what is fixed or changed
Added linear first-order constant-coefficient non-homogeneous system of ODEs solver.

#### Other comments

The SymPy's ode module will be able to solve a number of linear first-order constant coefficient non-homogenous systems of ODEs

For example:
```
>>> from sympy import *
>>> from sympy.abc import *
>>> f, g = symbols("f g", cls=Function)
>>> eqs = [Eq(f(t).diff(t), f(t) + 2*g(t) + t), Eq(g(t).diff(t), 7*g(t) + t**2)]
>>> dsolve(eqs)
[Eq(f(t), C1*exp(t) + C2*exp(7*t)/3 + 2*t**2/7 - 17*t/49 - 115/343), Eq(g(t), C2*exp(7*t) - t**2/7 - 2*t/49 - 2/343)]
```

#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your ease notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* solvers
  * In `dsolve` there is no a general solver that can solve systems of constant coefficient non-homogeneous first order ODEs of any size in terms of integrals.
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

mijo2 added 2 commits May 18, 2020 13:20
Updated the matching function to classify constant coefficient
non-homogeneous functions as type2 and dsolve uses the new solver.
Moved the test cases from test_ode.py to test_systems.py
@mijo2
Copy link
Contributor Author

mijo2 commented May 20, 2020

@oscarbenjamin @moorepants Except the final travis CI build and adding of few more test cases, the PR is almost completed. Please review and have a look. @namannimmo10 I hope this helps for your gsoc project.

mijo2 added 3 commits May 20, 2020 13:50
This commit may be squashed later.
Some of the matching function code was incorrectly removed and it
was rectified in this PR.
@mijo2
Copy link
Contributor Author

mijo2 commented May 21, 2020

@oscarbenjamin I don't know for what reason but the latex build is failing and I wasn't able to find out whats the problem.

@oscarbenjamin
Copy link
Contributor

Have you tried building the latex docs yourself?

This is the end of the build log:

LaTeX Warning: Citation `modules/solvers/ode:r694' on page 1202 undefined on in
put line 90359.
! Missing } inserted.
<inserted text> 
                }
l.90406 ... e^{- A t} b\, dt\right + C)\end{split}

Searching for that in the diff here I find:

.. math::
X = e^{A t} (\int e^{- A t} b\, dt\right + C)

I don't think that the \right should be there. The way to use \right is like \left( \cos{x} \right).

@mijo2
Copy link
Contributor Author

mijo2 commented May 21, 2020

Have you tried building the latex docs yourself?

I did but I wasn't able to figure the error out.

I don't think that the \right should be there. The way to use \right is like \left( \cos{x} \right).

I see, when I separately generated the equation, it worked and gave the right answer. Thanks for this help.


e1, e2, e3 = symbols("e1:4", real=True)
b1, b2, b3 = symbols("b1:4", real=True)
v1, v2, v3 = symbols("v1:4", cls=Function, real=True)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Rather than duplicating this code here and above I think that it would be better to make a function that returns the system and expected solution and then call that in each of the two test functions here.

eq1 = r1*c1*Derivative(x1(t), t) + x1(t) - x2(t) - r1*i
eq2 = r2*c1*Derivative(x1(t), t) + r2*c2*Derivative(x2(t), t) + x2(t) - r2*i
eq = [eq1, eq2]
sol = [Eq(x1(t), -2*C1*c2*r2*exp(-t/(2*c2*r2) - t/(2*c2*r1) - t/(2*c1*r1) - t*sqrt(c1**2*r1**2 +
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This can be simplified using the script.

Some of the test cases used repeated code, this was rectified in
this commit.
The big test cases were updated with smaller representations of their
solutions.
I = Function("I")
system = [Eq(V(t).diff(t), -1 / RC * V(t) + I(t) / C), Eq(I(t).diff(t), -R1 / L * I(t) - 1 / L * V(t) + Vs / L)]

return system, dsolve(system), [RC, t, C, Vs, L, R1, V0, I0, V, I]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be better to define the expected solution here and return system, funcs, sol and let the caller call dsolve. The way this is set up dsolve gets called twice because of the two tests that call this function.


eqs, dsolve_sol, symbols = _de_lorentz_solution()
m, q, t, e1, e2, e3, b1, b2, b3, v1, v2, v3 = symbols
x_1 = Integral(b1**3*b2*e2*q*exp(q*t*sqrt(-b1**2 - b2**2 - b3**2)/m)/(2*b1**3*b3*m +
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This long expression for x_1 can be simplified further by running printsol on it.

@@ -795,16 +1120,23 @@ def test_linear_3eq_order1_type4_skip():
Eq(diff(y(t), t), 2 * f * x(t) + (f + g) * y(t) - 2 * f * z(t)), Eq(diff(z(t), t), 5 * f * x(t) + f * y(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function shouldn't be called _skip if it's not being skipped.

Updated one of the big test cases so that we call the dsolve method
less number of times.
x_3 = exp(-t / (2 * c2 * r2) - t / (2 * c2 * r1) - t / (2 * c1 * r1) + t * sqrt(x_1) / (2 * c1 * c2 * r1 * r2))
x_4 = Integral(-i * r2 * exp(
t / (2 * c2 * r2) + t / (2 * c2 * r1) + t / (2 * c1 * r1) + t * sqrt(x_1) / (2 * c1 * c2 * r1 * r2)) / sqrt(
x_1), t)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some of these lines are too long. They wouldn't need to be so long though if you didn't put all the spaces in:

x_4 = Integral(-i*r2*exp(t/(2*c2*r2) + t/(2*c2*r1) + t/(2*c1*r1) + t*sqrt(x_1)/(2*c1*c2*r1*r2))/sqrt(x_1), t)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This example can be written as

_x1 = c1**2*r1**2 + 2*c1**2*r1*r2 + c1**2*r2**2 - 2*c1*c2*r1*r2 + 2*c1*c2*r2**2 + c2**2*r2**2
_x2 = Integral(i*r2*exp(-sqrt(_x1)*t/(2*c1*c2*r1*r2) + t/(2*c2*r2)
                            + t/(2*c2*r1) + t/(2*c1*r1))/sqrt(_x1), t)
_x3 = Integral(-i*r2*exp(sqrt(_x1)*t/(2*c1*c2*r1*r2) + t/(2*c2*r2)
                            + t/(2*c2*r1) + t/(2*c1*r1))/sqrt(_x1), t)
_x4 = exp(sqrt(_x1)*t/(2*c1*c2*r1*r2) - t/(2*c2*r2) - t/(2*c2*r1) - t/(2*c1*r1))
_x5 = exp(-sqrt(_x1)*t/(2*c1*c2*r1*r2) - t/(2*c2*r2) - t/(2*c2*r1) - t/(2*c1*r1))
sol = [
    Eq(x1(t),
        - 2*C1*_x5*c2*r2/(sqrt(_x1) + c1*r1 + c1*r2 - c2*r2)
        - 2*C2*_x4*c2*r2/(-sqrt(_x1) + c1*r1 + c1*r2 - c2*r2)
        - 2*_x2*_x4*c2*r2/(-sqrt(_x1) + c1*r1 + c1*r2 - c2*r2)
        - 2*_x3*_x5*c2*r2/(sqrt(_x1) + c1*r1 + c1*r2 - c2*r2)),
    Eq(x2(t), C1*_x5 + C2*_x4 + _x2*_x4 + _x3*_x5),
]

x_2 = 1/(2*b1**2*m*exp(q*t*sqrt(-b1**2 - b2**2 - b3**2)/m) + 2*b2**2*m*exp(q*t*sqrt(-b1**2 - b2**2 - b3**2)/m) +
2*b3**2*m*exp(q*t*sqrt(-b1**2 - b2**2 - b3**2)/m))
x_3 = exp(q*t*sqrt(-b1**2 - b2**2 - b3**2)/m)
x_4 = sqrt(b1**2 + b2**2 + b3**2)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This example can be written as

x1 = sqrt(b1**2 + b2**2 + b3**2)
x2 = exp(q*t*sqrt(-b1**2 - b2**2 - b3**2)/m)
x3 = 1/(b1**2*x1*x2 + b2**2*x1*x2)
x4 = 1/(2*b1**2*m*x2 + 2*b2**2*m*x2 + 2*b3**2*m*x2)
x5 = Integral(
        b1*b3*e1*q/(b1**2*m + b2**2*m + b3**2*m)
      + b2*b3*e2*q/(b1**2*m + b2**2*m + b3**2*m)
      + b3**2*e3*q/(b1**2*m + b2**2*m + b3**2*m), t)
x6 = 1/(2*b1**3*b3*m + 2*I*b1**2*b2*m*x1 + 2*b1*b2**2*b3*m
      + 2*b1*b3**3*m + 2*I*b2**3*m*x1 + 2*I*b2*b3**2*m*x1)
x7 = Integral(
        b1**2*e3*q*x4 - b1*b3*e1*q*x4 + I*b1*e2*q*x1*x4
      + b2**2*e3*q*x4 - b2*b3*e2*q*x4 - I*b2*e1*q*x1*x4, t)
x8 = Integral(
        b1**3*b2*e2*q*x2*x6 - b1**2*b2**2*e1*q*x2*x6 - b1**2*b3**2*e1*q*x2*x6
      - I*b1**2*b3*e2*q*x1*x2*x6 + b1**2*e3*q*x2/(2*b1**2*m + 2*b2**2*m + 2*b3**2*m)
      + b1*b2**3*e2*q*x2*x6 - b2**4*e1*q*x2*x6 - b2**2*b3**2*e1*q*x2*x6
      - I*b2**2*b3*e2*q*x1*x2*x6 + b2**2*e3*q*x2/(2*b1**2*m + 2*b2**2*m + 2*b3**2*m), t)
sol = [
    Eq(v1(t),
        C1*b1/b3
      - I*C2*b1**2*b2*x3 - C2*b1*b3*x1*x3 - I*C2*b2**3*x3 - I*C2*b2*b3**2*x3
      + I*C3*b1**2*b2*x2/(b1**2*x1 + b2**2*x1) - C3*b1*b3*x1*x2/(b1**2*x1 + b2**2*x1)
      + I*C3*b2**3*x2/(b1**2*x1 + b2**2*x1) + I*C3*b2*b3**2*x2/(b1**2*x1 + b2**2*x1)
      + I*b1**2*b2*x2*x7/(b1**2*x1 + b2**2*x1) - I*b1**2*b2*x3*x8
      - b1*b3*x1*x2*x7/(b1**2*x1 + b2**2*x1) - b1*b3*x1*x3*x8 + b1*x5/b3
      + I*b2**3*x2*x7/(b1**2*x1 + b2**2*x1) - I*b2**3*x3*x8
      + I*b2*b3**2*x2*x7/(b1**2*x1 + b2**2*x1) - I*b2*b3**2*x3*x8),
    Eq(v2(t),
        C1*b2/b3
      + C2*b1*sqrt(-b1**2 - b2**2 - b3**2)/(b1**2*x2 + b2**2*x2)
      - C2*b2*b3/(b1**2*x2 + b2**2*x2) - I*C3*b1*x1*x2/(b1**2 + b2**2)
      - C3*b2*b3*x2/(b1**2 + b2**2) - I*b1*x1*x2*x7/(b1**2 + b2**2)
      + b1*x8*sqrt(-b1**2 - b2**2 - b3**2)/(b1**2*x2 + b2**2*x2)
      - b2*b3*x2*x7/(b1**2 + b2**2) - b2*b3*x8/(b1**2*x2 + b2**2*x2) + b2*x5/b3),
    Eq(v3(t), C1 + C3*x2 + x2*x7 + x5 + (C2 + x8)*exp(-q*t*sqrt(-b1**2 - b2**2 - b3**2)/m)),
]

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok I will update this.

@oscarbenjamin
Copy link
Contributor

I've put the simplifying script in #19574 and improved it a bit.

@mijo2 mijo2 force-pushed the const_coeff_nonhomogeneous branch from b76beb1 to fac6758 Compare June 18, 2020 09:56
@mijo2
Copy link
Contributor Author

mijo2 commented Jun 18, 2020

@oscarbenjamin All checks have passed.

+ 2 * b1 * b3 ** 3 * m + 2 * I * b2 ** 3 * m * x1 + 2 * I * b2 * b3 ** 2 * m * x1)
x7 = Integral(
b1 ** 2 * e3 * q * x4 - b1 * b3 * e1 * q * x4 + I * b1 * e2 * q * x1 * x4
+ b2 ** 2 * e3 * q * x4 - b2 * b3 * e2 * q * x4 - I * b2 * e1 * q * x1 * x4, t)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand why you've added all of these spaces. Is that something that your editor does automatically? It's easier to read with the spacing from str and careful line breaks:

x7 = Integral(b1**2*e3*q*x4 - b1*b3*e1*q*x4 + I*b1*e2*q*x1*x4
            + b2**2*e3*q*x4 - b2*b3*e2*q*x4 - I*b2*e1*q*x1*x4, t)

Copy link
Contributor Author

@mijo2 mijo2 Jun 18, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Most probably it's the problem with my editor, I would get it fixed.

@oscarbenjamin
Copy link
Contributor

I would merge this but github won't let me. It hasn't updated since Travis passed the tests.

@oscarbenjamin
Copy link
Contributor

I've edited the release note and closed/reopened to restart the tests. Hopefully github will detect when the tests are complete on Travis (#19589)

@oscarbenjamin
Copy link
Contributor

I'm still unable to merge this even though the tests have passed (again)...

@oscarbenjamin
Copy link
Contributor

Apparently commenting updated the PR and made it mergeable!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
4 participants