Improved ODE Solver

[Manoj Kumar]

Hello Everyone,

I've been trying to solve bugs related to the ODE Solver .I would like my proposal to be in this area  since I've got a general idea of the source code. I found the following features missing in ode.py. I could possibly integrate these ideas into a single project called "Improving ODE solver in sympy" for my SOC proposal.

Most of my references are from Differential Equations with Applications abd Historical Notes from George F. Simmons

1 . Solving a set of homogenous equations with constant coefficients:
 
Right now sympy has no support for solving a system of differential equations. I would like to implement a method by which a two system diff. equation can be solved when there are constant coefficients. If a function of t is also present in the system of differential equations, then try using variation of parameters to solve the system.The actual session would be something like this. I thought of trying to implement it in dsolve, but wouldn't it make things a bit messy? So maybe use another function called dsolve_system(?)

>>> from sympy import *
>>> from sympy.abc import x, y, t
>>> f1 = Derivative(x, t)
>>> f2 = Derivative(y, t)
>>> eq1 = f1 - 3*x - (4*y)
>>> eq2 = f2 - x + y
>>> dsolve_system(eq1, eq2)
Output: [x = 2 * C1 * exp(t) + C2 * (1 + 2*t) * exp(t) , y = C1 * exp(t) + C2 * t * exp(t)]  
>>> eq1 = f1 + 3*sin(x) - (4*y)
>>> dsolve_system(eq1, eq2)
Output: ValueError: dsolve_system can work only when there are constant coefficients.

I have studied methods to implement this only when there are constant coefficients, if there are methods to do this when the coefficients are functions of x and y do let me know i'll read them up.

2.Power series solution of homogenous second order linear equations:
The user can input the point about which the power series solution can be found. If not the point can be assumed to be zero.

Implementing it in dsolve might be a pain, because the user would want the equation upto n terms and sometimes the nth term.
So, define a SeriesSolve class, evaluate a general term and return when a method called term is called (The simplest way I could think of, feel free to criticize me)
Something like this, maybe

>>> from sympy.solvers.ode import SeriesSolve
Lets look at various cases.
>>> f = Function('f')(x)
>>> d1 = Derivative(f, x)
>>> d2 = Derivative(f, x, x)
>>> eq = (1 - x**2)*d2 - (2*x)d1 + 6 * f
>>> series = SeriesSolve(eq)
>>> series.term(0)
Output: [A0 , A1]
>>> series.term(2)
Output: [-3 *A0 , - 2 * A1 / 3)
>>> series = SeriesSolve(eq, point = 2)
>>> series.term(0)
Output : [A0 , A1]
>>> series.term(2)
Output: [A0 , -2*A1/3]
>>> series.equation(n)
#Can give the first n terms of both the series.

And if the given point is singular check if it is regular, If regular attempt to find the Frobenius series solution.
>>> eq = (2*x**2)*d2 + x * (2*x + 1)*d1 - f
>>> series = SeriesSolve(eq)
>>> series.term(1)
Output: [(-2 / 5)*A0 , - A1]
>>> series.equation(n)
Same goes, first n terms of both the series

If the roots of the indicial equation are equal, return a ValueError. If the second series cannot be found, just return the first series.

If point is not regular, return ValueError
>>> eq = (2*x**3)*d2 + x * (2*x + 1)*d1 - f
Output: ValueError: Point is irregular singular, power series solution cannot be found.

And if the equation is of the form of a Gauss HyperGeometric one, apply the short-cut , transform it into a GHE , and get the general solution.

If possible try and implement initial value conditions. Since till now there seems to be no support for values like f'(0) I suppose the easiest way would be

>>> equation = SeriesSolve(eq, point = x, (point1 , value of function) , (point2, value of derivative of function))

* This would be a bit tough in the case where there is a three recursion formula, where I'm not sure but I shall verify this as soon as possible.

* In cases like Frobenius Series Solution, sometimes there might be only one series, then simply ignore the second condition?

3.Implementing Laplace transform to solve equations of second degree with initial conditions.

Laplace transforms could be used in a simpler way to solve the above mentioned equations.

This could be again implemented in dsolve, but since initial conditions aren't specified, this could again be a big pain.

>>> from sympy.solvers.ode import laplace_eq
>>> laplace_eq(d2 - 4*d1 + 4*f , (0 , 0) , (0, 3))
Output:
3 * x * exp(2*x)

This could be extended to nth order, but supplying initial conditions and finding laplace inverse would then become difficult

These are the ideas I had done as part of my coursework last semester. Please do tell me , if they are trivial(which I hope not) , places where I could improve, maybe new methods that I could try and read up to fit in a period of ten weeks , so that the possibilty of me getting selected could increase.

Regards,
Manoj Kumar,
Mech Undergrad.
BPGC
Blog

[Stefan Krastanov]

Hi

Concerning only

>1 . Solving a set of homogenous equations with constant coefficients:

check out https://github.com/sympy/sympy/pull/1322

There you can find a lot a discussion about how it should be done and
code that solves not only inhomogeneous linear systems but also some
nonlinear ones. It breaks on degenerate cases, which can not yet be
done, because the matrix module does not support certain
generalizations of eigenvalue decomposition.

It is quite problematic to write tests for this module because of the
mess with integration constants.

[Stefan Krastanov]

If you are interested in the project you will have to solve to
substantial issues (regardless of whether you use the code that I
suggested):

- creating an actual "integration constant" object (and making all
the testing code to work with it)
- adding new decomposition routines to the matrix module

[someone]

> - adding new decomposition routines to the matrix module

Actually implementing Jordan decompositions and matrix exponentials I think.

[Manoj Kumar]

Thanks Stefan,

I've solved just homogeneous system of equations, and do the problems you speak of (seems a bit higher level math to me) exist in homogeneous equations also?

Could you suggest a few docs, or books that I should read , so that I get my algorithm clear before looking into your code?

 

On Sun, Feb 24, 2013 at 6:59 PM, someone <some...@bluewin.ch> wrote:

>  - adding new decomposition routines to the matrix module

Actually implementing Jordan decompositions and matrix exponentials I think.

--
You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/C4btttnGJss/unsubscribe?hl=en.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.

--

[Manoj Kumar]

And since you have already done parts of it, what is the scope for me to do it?

[Stefan Krastanov]

concerning work to be done: I have created only a skeleton
ode_systems.py solver. While I like the design, it is incomplete (it
raises NotImplementedError) and not at all integrated in the rest of
sympy. And all the rest of the projects suggested by you are still
possible.

The parts that are incomplete will result in only a few lines of code
_if_ we have jordan decomposition.

If you know what eigenvalue decomposition is and why it is necessary
for solving systems of linear differential equations, then the rest is
a straightforward generalization for solving systems that can not be
decomposed in this way. Check the code that I have written in
ode_systems.py (it is commented) and wikipedia, and work out a few
simple examples by hand (check the failing tests and try them step by
step with sage/mathematica/wolfram alpha)

If you wish you can cannibalize my code and provide a pull request in
whatever style you want.

[Manoj Kumar]

Thanks, Sorry for bugging but it would be nice if other mentors could also give their review or feedback on my ideas, especially the second one.

[Manoj Kumar]

Regarding my second idea,

I've gone through the sympy documentation and found out that at present there is also no support, to identify whether a given function is regular or irregular singular at a point or an interval. And I've also found about that being a GSOC idea in itself. i.e Analysis of function in a given interval.

Since it is impossible to get a series solution, without actual analysis of the respective functions (Coefficients of y'' , y' and y)at the point, I guess that would also come under my second project idea. It is also given in the idea list that it is a bit vague and needs discussion in the group.So maybe I could stick on only to the second idea, if it does seem enough for a period of 12 weeks.

Awaiting your feedback.

[Aaron Meurer]

As Stefan noted, there is already some work on this. In general, you
can solve a system of n equations with n variables with constant
coefficients in full generality. When this sort of thing happens, it
is always better to implement the full solution rather than just a
special case. For example, you may have learned how to apply
variation of parameters for second order ODEs, but in fact it can be
applied to nth order ODEs, and this is the method implemented in
SymPy.

Actually, we do support that:

In [10]: print f(x).diff(x).subs(x, 0)
Subs(Derivative(f(x), x), (x,), (0,))

In [11]: print f(x).diff(x).subs(x, a)
Derivative(f(a), a)

When I originally wrote the module, there wasn't support for the Subs
object, so writing f'(0) was indeed impossible, which is why initial
condition support was not added, but now that there is, it should not
be too hard to add it.

See issue https://code.google.com/p/sympy/issues/detail?id=1621 and
also https://code.google.com/p/sympy/issues/detail?id=1620 for much
discussion on this.

>
>>>> equation = SeriesSolve(eq, point = x, (point1 , value of function) ,
>>>> (point2, value of derivative of function))
>
> * This would be a bit tough in the case where there is a three recursion
> formula, where I'm not sure but I shall verify this as soon as possible.
>
> * In cases like Frobenius Series Solution, sometimes there might be only one
> series, then simply ignore the second condition?
>
> 3.Implementing Laplace transform to solve equations of second degree with
> initial conditions.
>
> Laplace transforms could be used in a simpler way to solve the above
> mentioned equations.
>
> This could be again implemented in dsolve, but since initial conditions
> aren't specified, this could again be a big pain.
>
>>>> from sympy.solvers.ode import laplace_eq
>>>> laplace_eq(d2 - 4*d1 + 4*f , (0 , 0) , (0, 3))
> Output:
> 3 * x * exp(2*x)
>
> This could be extended to nth order, but supplying initial conditions and
> finding laplace inverse would then become difficult
>
> These are the ideas I had done as part of my coursework last semester.
> Please do tell me , if they are trivial(which I hope not) , places where I
> could improve, maybe new methods that I could try and read up to fit in a
> period of ten weeks , so that the possibilty of me getting selected could
> increase.

In some sense, some of these are trivial, because a lot of the hard
work has already been done. For example, there is already support for
handling laplace transforms with the function laplace_transform(), so
one would just need to write up the algorithm using this function, and
iron out any issues that arise. Regarding the series stuff, this
would be more work. As I recall, in general, you get a recurrence
relation on the coefficients. I don't know just how powerful our
rsolve() is, but it seems to do alright, so that in many cases at
least you can generate an exact summation for the solution (as opposed
to an expansion up to a given but finite order). Your idea for an
object to represent a series solution is good thinking. We will
probably need something like that in the general case.

You can always find more ODE solving method to implement. Loop up all
the methods that Maple uses in their docs, for example. In general,
for such a project, you'll end up learning a lot of things, so don't
just try to limit yourself to what you already know.

Aaron Meurer

>
> Regards,
> Manoj Kumar,
> Mech Undergrad.
> BPGC
> Blog
>

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an

[Manoj Kumar]

Aaron: Thanks for your feedback.

I've seen the maple docs, and found out that they use lie groups and symmetry methods also to solve ODE's. It is not implemented in sympy. I've just started reading up on it. Do you think it would be a good idea to include it in my SOC proposal?

You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/C4btttnGJss/unsubscribe?hl=en.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.

--

[Stefan Krastanov]

> I've seen the maple docs, and found out that they use lie groups and
> symmetry methods also to solve ODE's. It is not implemented in sympy. I've
> just started reading up on it. Do you think it would be a good idea to
> include it in my SOC proposal?

If done well, such additional solvers (or rather sub-solvers) may even
be sufficient for a project on their own. It all depends on how you
structure the application. More suggested functionality is certainly
better if you can show that you can deliver on you promises.

One thing to keep in mind is that the more advance the math you will
be working with, the more important it gets for you to actually prove
that you know the math. This can be done in may ways, in the
application or in a preliminary pull request.

[Manoj Kumar]

Hello again.

Sorry to bring this post up again, but I thought it would be good to keep all my ideas on a single thread.

I've been reading about lie groups and their application in solving differential equations and I believe I've got my concepts right about solving a linear first order differential equation. Maple, (according to the papers cited in the SoC ideas pages), has two separate papers dedicated to solving differential equations of the first and second order using lie groups. I would like to like to code a clone for solving first order differential equations using lie groups in sympy for my SoC proposal.

Theory:
There are many differential equations that cannot be solved using the present dsolve in ode.py. Suppose there is a differential equation of the form
(dr / ds) = f(r), since the slope is independent of s, moving in the s direction would merge one solution into another. The principal aim of using lie groups is to reduce any given differential equation to the above mentioned form.

The coordinates r and s are called canonical coordinates because any lie group transformation which leaves the differential equation invariant would transform it into the form (r, s + alpha) where r remains invariant.If we are able to find such coordinates, then the equation becomes separable.

The problem is finding these canonical co-ordinates. This can be broken into four steps(This is for first order differential equations)

1. Find η and ξ , the infinitesimals using this partial diffential equation, which
are basically the differentials of the transformed co-ordinates with respect to
the transformation parameter.

ηx − ξy h2 + (ηy − ξx )h − (ξhx + ηhy ) = 0

This should be the toughest part of the project. As you can see solving a linear first order ODE, leads to solving a much more complicated PDE. However in this PDE, some assumptions could be made for the infinitesimals.

Maple has 6 separate algorithms just for finding the infinitesimals as described in the paper. Once this is done.

2. rx ξ + ry η = 0,
   sx ξ + sy η = 1,

Solving these two equations would give, r and s , this can be done by. dx / ξ = dy / η = r(x,y) and dx / ξ = ds

3. After this substituting r and s  in
ds / dr = sx + h(x,y)sy / rx + h(x, y)ry and simplifying in terms of r and s would give a separable form.

4. And after solving the differential equation in terms of r and s, substituting back r and s in terms of x and y.

I haven't focussed on the implementation part yet, but I think it should be structured similar to dsolve, or should it be implemented in dsolve.

I've focussed on only first order differential equations here, even Maple also has support only upto the second order.I feel this project would be challenging for me and also a addition to the ODE solver of sympy.

@Aaron, Stefan, and Smichr Please do give your feedback.

[Chris Smith]

> @Aaron, Stefan, and Smichr Please do give your feedback.

I defer to those who work and study these things.

[Aaron Meurer]

> ηx - ξy h2 + (ηy - ξx )h - (ξhx + ηhy ) = 0

>
> This should be the toughest part of the project. As you can see solving a
> linear first order ODE, leads to solving a much more complicated PDE.
> However in this PDE, some assumptions could be made for the infinitesimals.
>
> Maple has 6 separate algorithms just for finding the infinitesimals as
> described in the paper. Once this is done.
>
> 2. rx ξ + ry η = 0,
> sx ξ + sy η = 1,
>
> Solving these two equations would give, r and s , this can be done by. dx /
> ξ = dy / η = r(x,y) and dx / ξ = ds
>
> 3. After this substituting r and s in
> ds / dr = sx + h(x,y)sy / rx + h(x, y)ry and simplifying in terms of r and s
> would give a separable form.
>
> 4. And after solving the differential equation in terms of r and s,
> substituting back r and s in terms of x and y.
>
> I haven't focussed on the implementation part yet, but I think it should be
> structured similar to dsolve, or should it be implemented in dsolve.
>
> I've focussed on only first order differential equations here, even Maple
> also has support only upto the second order.I feel this project would be
> challenging for me and also a addition to the ODE solver of sympy.
>
> @Aaron, Stefan, and Smichr Please do give your feedback.

I haven't formally studied Lie groups and their applications to
solving ODEs, but from my understanding, they are indeed a powerful
method for doing so. This would definitely be a welcome project for
GSoC.

Aaron Meurer

[Manoj Kumar]

Thanks for your reply. So do I start formally working on my proposal (Solving first order ODE's using Lie Groups), that is creating a wiki page for it. And if I do, who is most likely to be my mentor? . Because I've seen quite a few good proposals not accepted last time, because of a need of a mentor.

You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/C4btttnGJss/unsubscribe?hl=en.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.

[Aaron Meurer]

Don't worry about not being accepted for this reason. If we decide we
want the project, we'll find someone to mentor it. It's true that we
are limited by how many proposals we accept mostly by how many mentors
we have, but the main thing that should concern you is making the
proposal high quality. And anyway, this proposal in particular I
think almost anyone could mentor should it come down to it. That is
really only an issue for proposals that involve very heavy math (like
some project in the polys), or some physics (like the quantum or
mechanics module), where we need to make sure the mentor understands
everything.

Aaron Meurer

On Fri, Mar 8, 2013 at 12:35 AM, Manoj Kumar

[someone]

> (Solving first order ODE's using Lie Groups)

As far as I read some papers about this, the Lie methods
become most powerful in the 2nd (and maybe higher) order
case. So I suggest that you consider also 2nd order
problems. No question any improvement on ODE solvers
is highly desirable.

[Manoj Kumar]

@someone: Thanks a lot for your suggestion.

However, I believe even for first order differential equations, there is a lot of work to be done. For solving the partial differential equation there are a set of six algorithms, each having certain sub-algorithms , which I guess in itself would take 1 - 1.5 weeks each to implement. I could always add support to 2nd order ODE's as well if there is sufficient time left.

Also, I would like to know if the following suggestions are right.

1. It should be implemented as a hint in dsolve itself, which would be the best possible way

2. I should make a separate module to find the infinitesimals , that is If I give the differential equation, it should return the infinitesimals. This would be good, if someone wants to add support to 2nd order ODE's and PDE's later on.

[Mary Clark]

I think that a project like this could tie in quite well with a separate Lie group/lie algebra module, especially if we have something that would store a lot of information about some of the most common Lie groups (SL(2,C), SU(2), SO(3), etc). 

[Manoj Kumar]

Hello,

Pardon me for my ignorance, but my knowledge of lie groups is limited to the solving of differential equations. However it seem a really good idea.

And I don't think that actual finding of lie groups is used in differential equations, it is just the application that is used in finding infinitesimals, as mentioned in my post which are used for making transformations.

On Fri, Mar 8, 2013 at 8:27 PM, Mary Clark <mary.sp...@gmail.com> wrote:

I think that a project like this could tie in quite well with a separate Lie group/lie algebra module, especially if we have something that would store a lot of information about some of the most common Lie groups (SL(2,C), SU(2), SO(3), etc). 

--

[Aaron Meurer]

On Mar 8, 2013, at 4:51 AM, Manoj Kumar <manojkumar...@gmail.com> wrote:

@someone: Thanks a lot for your suggestion.

However, I believe even for first order differential equations, there is a lot of work to be done. For solving the partial differential equation there are a set of six algorithms, each having certain sub-algorithms , which I guess in itself would take 1 - 1.5 weeks each to implement. I could always add support to 2nd order ODE's as well if there is sufficient time left.

Also, I would like to know if the following suggestions are right.

1. It should be implemented as a hint in dsolve itself, which would be the best possible way

Yes, the user level interface should be dsolve. Most users don't care how it is done. They just want to pass the ode to dsolve and get a solution. 

2. I should make a separate module to find the infinitesimals , that is If I give the differential equation, it should return the infinitesimals. This would be good, if someone wants to add support to 2nd order ODE's and PDE's later on.

I don't fully understand what it entails yet, so it's hard to say. However, it's easy to move code around, especially internal code, so I wouldn't worry about this yet. 

Aaron Meurer

 

So I suggest that you consider also 2nd order

problems.

--