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

Exact ODE solver class added #18887

Open
wants to merge 12 commits into
base: master
Choose a base branch
from
Open

Exact ODE solver class added #18887

wants to merge 12 commits into from

Conversation

tnzl
Copy link

@tnzl tnzl commented Mar 17, 2020

References to other Issues or PRs

Refactoring efforts in accordance with #18348.

Brief description of what is fixed or changed

Exact ODE solver has been added with conditions to convert given ODE to exact where possible.
ode_classify() in ode.py has been updated to use the FirstExact solver class.
_handle_Integral() in ode.py has been refactored to remove global y hack.

Other comments

Release Notes

  • solvers

    • Added FirstExact solver class for the exact ordinary differential equation.
  • functions

    • Refactored ode_classify() and _handle_Integral() in ode.py

@sympy-bot
Copy link

sympy-bot commented Mar 17, 2020

Hi, I am the SymPy bot (v158). 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:

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

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
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Refactoring efforts in accordance with #18348.

#### Brief description of what is fixed or changed
Exact ODE solver has been added with conditions to convert given ODE to exact where possible.
ode_classify() in ode.py has been updated to use the FirstExact solver class.
_handle_Integral() in ode.py has been refactored to remove global y hack.

#### Other comments


#### 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 release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* solvers
  * Added FirstExact solver class for the exact ordinary differential equation.

* functions
  * Refactored ode_classify() and _handle_Integral() in ode.py
<!-- END RELEASE NOTES -->

y=Dummy('y')

self.q=collect(eq,f.diff(x)).coeff(f.diff(x))
self.p=eq-self.q*f.diff(x)
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 use the pattern matching approach like FirstPatternSolver

Copy link
Author

Choose a reason for hiding this comment

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

Thanks for the suggestion, I will change it to inherit SinglePatternODESolver.

def _matches(self) -> bool:
if self.ode_problem.order!=1:
self._matched=False
return False
Copy link
Contributor

Choose a reason for hiding this comment

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

It isn't necessary to set _matched in this function

Copy link
Author

Choose a reason for hiding this comment

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

Yeah! matches() sets _matched. Updating it.

hint="1st_exact"
has_integral=True
p=None
q=None
Copy link
Contributor

Choose a reason for hiding this comment

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

There should be spaces around =

class FirstExact(SingleODESolver):
r"""Solves 1st order exact ordinary differential equations.
A 1st order differential equation is called exact if it is the total
differential of a function. That is, the differential equation:
Copy link
Contributor

Choose a reason for hiding this comment

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

There already is a docstring for this method in ode_1st_exact. The ode_1st_exact function should be removed but the docstring can be moved here

try:
sol=sol.doit()
gen_sol=Eq(sol.subs(y,f),C1)
except:
Copy link
Contributor

Choose a reason for hiding this comment

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

Bare except should not be used anywhere in the codebase

Copy link
Author

Choose a reason for hiding this comment

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

Removing try except block as doit will not be called here.


try:
sol=sol.doit()
gen_sol=Eq(sol.subs(y,f),C1)
Copy link
Contributor

Choose a reason for hiding this comment

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

We don't need to call doit here as it will be called by dsolve later.

Copy link
Author

Choose a reason for hiding this comment

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

Updating it. As there is no other way to pass y to _handle_Integral, global y hack will still be in use.

try:
if simplify(self.p)!=0:
m=self.p.subs(f,y)
n=self.q.subs(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.

There should be spaces after comma (see PEP8)

Copy link
Author

Choose a reason for hiding this comment

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

Thanks for the suggestion, I will format to PEP8 style.

@oscarbenjamin
Copy link
Contributor

I'm not sure the global y is needed at all because we can integrate with functions as integration variables:

In [5]: Integral(f(x), (f(x), 0, 1))                                                                                              
Out[5]: 
1             
⌠             
⎮ f(x) d(f(x))
⌡             
0             

In [6]: Integral(f(x), (f(x), 0, 1)).doit()                                                                                       
Out[6]: 1/2

Maybe it's just not necessary to substitute f(x) for y in the first place. Can the hack be removed like that?

def _wilds(self, f, x, order):
P = Wild('P')
Q = Wild('Q', exclude=[f(x).diff(x)])
return P, Q
Copy link
Contributor

Choose a reason for hiding this comment

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

Should P also exclude f(x).diff(x)?

What happens with something like this?

1/f(x).diff(x) + f(x).diff(x)

Copy link
Author

Choose a reason for hiding this comment

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

Pattern match fails for equations which have Derivative in more than one term if P excludes f(x).diff(x), like eq=x*f(x).diff(x)+(2/f(x))*f(x).diff(x)+f(x).

Copy link
Author

Choose a reason for hiding this comment

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

Should P also exclude f(x).diff(x)?

What happens with something like this?

1/f(x).diff(x) + f(x).diff(x)

At present _matches checks equations of degree 1 only.

Copy link
Contributor

Choose a reason for hiding this comment

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

Maybe we should collect on derivatives here:

eq = eq.collect(f(x), func = cancel)

if fx.diff(x) in self._wilds_match[P].atoms(Derivative):
eq = collect(eq,fx.diff(x))
pattern = self._equation(fx, x, 1)
self._wilds_match = eq.match(pattern)
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 what this is doing but it doesn't seem right

Copy link
Contributor

Choose a reason for hiding this comment

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

The _wilds_match variable should be set by PatternSolver.match.

Copy link
Author

Choose a reason for hiding this comment

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

For the case where Derivative is present in more then one terms, eq needs to be preprocessed before matching. This condition is required to process and verify that P is the function of x , f(x) only and independent of any derivative.

Copy link
Contributor

Choose a reason for hiding this comment

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

I think that should be handled by excluding f(x).diff(x) in the pattern

Copy link
Author

Choose a reason for hiding this comment

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

The reason for not excluding f(x).diff(x) in Wild P and the need for this condition here are the same. Lets take an example of exact equation:

>>> P1 = Wild('P1',exclude=[f(x).diff(x)])
>>> Q1 = Wild('Q1',exclude=[f(x).diff(x)])
>>> 
>>> P2 = Wild('P2')
>>> Q2 = Wild('Q2',exclude=[f(x).diff(x)])
>>> 
>>> exprn1= P1 + Q1*f(x).diff(x)
>>> exprn2= P2 + Q2*f(x).diff(x)
>>> 
>>> eq=x*f(x).diff(x)+(2/f(x))*f(x).diff(x)+f(x)
>>> eq = eq.collect(f(x), func = cancel)
>>> 
>>> print(eq.match(exprn1))
None
>>> eq.match(exprn2)
{Q2_: 2/f(x), P2_: x*Derivative(f(x), x) + f(x)}

Here, match with exprn1 returns None. Hence, excluding f(x).diff(x) in Wild P will fail the matches function for for the case where f(x).diff(x) is present in more then one terms of eq.
After we match if f(x).diff(x) is present in P, we update it in above if block in _verify function.

Copy link
Author

@tnzl tnzl Mar 21, 2020

Choose a reason for hiding this comment

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

Every PatternODESolver will have different preprocess for matching. We need to define a _preprocess_eq function in the base class which preprocesses equation before it is matched in _matches function.

Copy link
Contributor

Choose a reason for hiding this comment

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

We need to collect on the derivatives of f as well:

In [7]: eq.collect(f(x).diff(x))                                                                                                  
Out[7]:2  ⎞ d              
⎜x + ────⎟──(f(x)) + f(x)
⎝    f(x)⎠ dx             

In [8]: eq.collect(f(x).diff(x)).match(exprn1)                                                                                    
Out[8]:2  ⎫
⎨P₁: f(x), Q₁: x + ────⎬
⎩                  f(x)⎭

Copy link
Author

Choose a reason for hiding this comment

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

Absolutely. We need to collect on the derivatives of f, but this is not required by other solvers. Therefore, putting eq.collect(f(x).diff(x)) in _matches of SinglePatternODESolver will not be a wise option.
We need _preprocess_eq in SinglePatternODESolver:

#Subclass should define here all preprocess on the equation required before pattern matching.
def _preprocess_eq(self, eq):
   prep_eq=eq
    return prep_eq

eq = _preprocess_eq(eq) will be called before eq.match(pattern) here:

self._wilds_match = match = eq.match(pattern)

After _preprocess_eq is defined in SinglePatternODESolver, we would implement it in FirstExact as:

def _preprocess_eq(self, eq):
    prep_eq = eq.collect(f(x).diff(x))
    return prep_eq

This will solve the problem of specific preprocessing required by different solvers before pattern matching

Copy link
Contributor

Choose a reason for hiding this comment

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

We need to minimise solver-specific logic. If collecting on f(x).diff(x) is generally useful then it should be done in the SingleODEPatternSolver class. Is there any reason not to put it there?

Normally collecting on f(x) also collects on derivatives e.g.:

In [60]: eq = f(x).diff(x) + x*f(x).diff(x) + f(x) + x*f(x)                                                                       

In [61]: eq                                                                                                                       
Out[61]: 
           d                 d       
xf(x) + x──(f(x)) + f(x) + ──(f(x))
           dx                dx      

In [62]: eq.collect(f(x))                                                                                                         
Out[62]: 
                       d       
(x + 1)f(x) + (x + 1)──(f(x))
                       dx 

The reason that doesn't work in your example is because you have mixed products involving both f(x) and f(x).diff(x):

In [63]: eq = f(x).diff(x) + (x/f(x))*f(x).diff(x) + f(x) + x*f(x)                                                                

In [64]: eq                                                                                                                       
Out[64]: 
           d                         
         x──(f(x))                  
           dx                d       
xf(x) + ────────── + f(x) + ──(f(x))
            f(x)             dx      

In [65]: eq.collect(f(x))                                                                                                         
Out[65]: 
  d                                 
x──(f(x))                          
  dx                        d       
────────── + (x + 1)f(x) + ──(f(x))
   f(x)                     dx      

In [66]: eq.collect(f(x)).collect(f(x).diff(x))                                                                                   
Out[66]: 
               ⎛ x      ⎞ d       
(x + 1)f(x) + ⎜──── + 1──(f(x))
               ⎝f(x)    ⎠ dx 

I think that the right way to do this is like this in the first place:

In [68]: eq.collect([f(x).diff(x), f(x)])                                                                                         
Out[68]: 
               ⎛ x      ⎞ d       
(x + 1)f(x) + ⎜──── + 1──(f(x))
               ⎝f(x)    ⎠ dx 

For ODE pattern matching collecting on the derivative is more important than collecting on f(x).

Copy link
Author

Choose a reason for hiding this comment

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

eq.collect([f(x).diff(x)], func = cancel) solves the problem for preprocess for all solvers including a second order Liouville solver for which I will open a PR separately.

factor = exp(Integral(factor).doit())
m*= factor
n*= factor
self._wilds_match[P] = m.subs(y,fx)
Copy link
Contributor

Choose a reason for hiding this comment

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

There should be spaces after commas and on both sides of something like *=. See PEP 8.

@tnzl
Copy link
Author

tnzl commented Mar 20, 2020

I'm not sure the global y is needed at all because we can integrate with functions as integration variables:

In [5]: Integral(f(x), (f(x), 0, 1))                                                                                              
Out[5]: 
1             
⌠             
⎮ f(x) d(f(x))
⌡             
0             

In [6]: Integral(f(x), (f(x), 0, 1)).doit()                                                                                       
Out[6]: 1/2

Maybe it's just not necessary to substitute f(x) for y in the first place. Can the hack be removed like that?

The general solution contains integrals that require f(x) and x to be treated as independent of each other.
This is the reason for substituting f(x) with y till the step where integration is done. At present there is no other way to pass y to _handle_Integral.

@oscarbenjamin
Copy link
Contributor

oscarbenjamin commented Mar 20, 2020

At present there is no other way to pass y to _handle_Integral.

Can the general solution be returned as a Subs?

e.g.:

In [45]: Subs(Integral(x*y, y), y, f(x))                                                                                          
Out[45]: 
⎛⌠       ⎞│      
⎜⎮ xy dy⎟│      
⎝⌡       ⎠│y=f(x)

In [46]: Subs(Integral(x*y, y), y, f(x)).doit()                                                                                   
Out[46]: 
   2   
xf (x)
───────
   2  

@tnzl
Copy link
Author

tnzl commented Mar 21, 2020

At present there is no other way to pass y to _handle_Integral.

Can the general solution be returned as a Subs?

e.g.:

In [45]: Subs(Integral(x*y, y), y, f(x))                                                                                          
Out[45]: 
⎛⌠       ⎞│      
⎜⎮ xy dy⎟│      
⎝⌡       ⎠│y=f(x)

In [46]: Subs(Integral(x*y, y), y, f(x)).doit()                                                                                   
Out[46]: 
   2   
xf (x)
───────
   2  

The general solution must be returned in integral form. Integral would be handled by dsolve later.

@oscarbenjamin
Copy link
Contributor

The general solution must be returned in integral form. Integral would be handled by dsolve later.

Can it be returned as Subs(Integral)?

@tnzl
Copy link
Author

tnzl commented Mar 21, 2020

Yes!!! This works perfectly.

collected on f(x).diff(x) before matching pattern

Updated wild variable P of FirstExact

Removed unused imports and variables
if simplify(m) != 0:
m = m.subs(fx,y)
n = n.subs(fx,y)
numerator= simplify(m.diff(y)-n.diff(x))
Copy link
Contributor

Choose a reason for hiding this comment

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

I still think that all of this code can be simplified somehow. Are all these calls to simplify necessary? In general we try to avoid calling simplify in library code because it can be slow. It is better to call specific simplify functions like factor, cancel etc. if we know what is needed.

Copy link
Contributor

Choose a reason for hiding this comment

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

Maybe we can just use something like if m.is_zero rather than if simplify(m) != 0.

Copy link
Author

Choose a reason for hiding this comment

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

simplify is called only before critical comparisons for converting to exact ODE, therefore removing doesn't seem possible. is_zero would not work here because it doesn't simplify expression before returning. I think this is the best way to simplify at present.

Copy link
Contributor

Choose a reason for hiding this comment

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

I would rather not call simplify during the matching code so it might be better just to be inconclusive in some cases.

Can you give an example of what m and n might be here so that simplify is needed?

Copy link
Author

@tnzl tnzl Mar 28, 2020

Choose a reason for hiding this comment

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

is_zero won't work here because it doesn't simplify the expression before return. Here, simplify is called only before critical conditions to convert to exact ODE therefore simplify can't be removed.
At present I think simplify is the best way to simplify expressions here.

Copy link
Contributor

Choose a reason for hiding this comment

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

Can you give the example of what m and n might be? There is almost always a better option then using generic simplify.

except NotImplementedError:
# Differentiating the coefficients might fail because of things
# like f(2*x).diff(x). See issue 4624 and issue 4719.
pass
Copy link
Contributor

Choose a reason for hiding this comment

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

This try block is too broad. There is too much code in between the try and the except. We should only catch exceptions around the narrowest block of code possible (ideally a single function call). Where does the NotImplemented error come from? Maybe the try/except can just be removed.

Copy link
Author

Choose a reason for hiding this comment

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

Are the tests failing because of my master branch being out of sync with sympy's master?

@oscarbenjamin
Copy link
Contributor

The tests failed because of #18978

@tnzl
Copy link
Author

tnzl commented Mar 29, 2020

@oscarbenjamin can you please guide me on why the tests are failing? All tests seems to have passed.

@oscarbenjamin
Copy link
Contributor

I think it was just a network problem. I've restarted the tests that stalled.

@codecov
Copy link

codecov bot commented Mar 29, 2020

Codecov Report

Merging #18887 into master will increase coverage by 0.017%.
The diff coverage is 96.551%.

@@              Coverage Diff              @@
##            master    #18887       +/-   ##
=============================================
+ Coverage   75.777%   75.795%   +0.017%     
=============================================
  Files          647       647               
  Lines       168623    168636       +13     
  Branches     39734     39732        -2     
=============================================
+ Hits        127778    127818       +40     
+ Misses       35280     35260       -20     
+ Partials      5565      5558        -7     

@tnzl
Copy link
Author

tnzl commented Mar 30, 2020

Thanks. Are there any other suggestions that you want me to implement here?

@oscarbenjamin
Copy link
Contributor

You haven't addressed the comments above

@oscarbenjamin
Copy link
Contributor

@tnzl what is the status of this?

@czgdp1807
Copy link
Member

@tnzl Are you still working on this?

@czgdp1807
Copy link
Member

@oscarbenjamin Can we close this or can it be taken over by some other contributor in future?

@oscarbenjamin
Copy link
Contributor

This is still needed so it would be good if someone could take this over.

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

Successfully merging this pull request may close these issues.

None yet

4 participants