[sympy] r3292 committed - Fixed Formatting
codesite...@google.com
Author: asmeurer
Date: Tue Oct 20 19:11:04 2009
Log: Fixed Formatting
http://code.google.com/p/sympy/source/detail?r=3292
Modified:
/wiki/ODEsModuleReport.wiki
=======================================
--- /wiki/ODEsModuleReport.wiki Sat Oct 10 09:51:47 2009
+++ /wiki/ODEsModuleReport.wiki Tue Oct 20 19:11:04 2009
@@ -1,9 +1,9 @@
#summary Review of the 2009 Summer of Code project
#labels GSoC
-*This page is still a work in progress. See my blog post for the full
report.*
-
-This review comes from my
[http://asmeurersympy.wordpress.com/2009/09/07/google-summer-of-code-2009-wrap-up/
blog]. See my blog for more details on my 2009 Google Summer of Code
experiences. Also, the blog post has some advice for anyone who may want
to apply for GSoC in the future.
+*See my
[http://asmeurersympy.wordpress.com/2009/09/07/google-summer-of-code-2009-wrap-up/
blog post] for the full report.*
+
+See my [http://asmeurersympy.wordpress.com/ blog] for more details on my
2009 Google Summer of Code experiences. Also, the blog post has some
advice for anyone who may want to apply for GSoC in the future.
=About Me=
My name is Aaron Meurer, and I am a second year math and computer science
major at New Mexico Tech in Socorro, New Mexico, USA.
@@ -52,30 +52,30 @@
==The GSoC Period==
For the start of the period, I followed my timeline. I implemented 1st
order ODEs with homogeneous coefficients and 1st order exact ODEs. These
were both pretty simple to do, as I expected.
-The next thing I wanted to do was separable. My goal at this point was to
get every relevant exercise from my textbook to work with my solvers. One
of the exercises from my <a
href="http://books.google.com.np/books?id=29utVed7QMIC&lpg=PA24&ots=uxLSUKt_3P&hl=en&pg=PA56#v=onepage&q=&f=false">book</a>
(Pg. 56, No. 21) was `dy/dx=exp**(x + y)`. I soon discovered that it was
impossible for SymPy to separate `e**(x + y)` into `exp(x)*exp(y)`, because
the second would be automatically combined in the core. I also discovered
that <code>expand()</code>, which should have been the function to split
that, expanded using all possible methods indiscriminately. Part of my
`separatevars()` function that I was writing to separate variables in
expressions would be to split things like `x + y*x` into `x*(1 + y)` and
`2*x*y + x**2 + y**2` into `(x + y)**2$, but `expand()`
+The next thing I wanted to do was separable. My goal at this point was to
get every relevant exercise from my textbook to work with my solvers. One
of the exercises from my <a
href="http://books.google.com.np/books?id=29utVed7QMIC&lpg=PA24&ots=uxLSUKt_3P&hl=en&pg=PA56#v=onepage&q=&f=false">book</a>
(Pg. 56, No. 21) was `dy/dx=exp**(x + y)`. I soon discovered that it was
impossible for SymPy to separate `e**(x + y)` into `exp(x)*exp(y)`, because
the second would be automatically combined in the core. I also discovered
that `expand()`, which should have been the function to split that,
expanded using all possible methods indiscriminately. Part of my
`separatevars()` function that I was writing to separate variables in
expressions would be to split things like `x + y*x` into `x*(1 + y)` and
`2*x*y + x**2 + y**2` into `(x + y)**2`, but `expand()`
as it was currently written would expand those.
-So I spent a few weeks hacking on the core to make it not auto-combine
exponents. I came up with a rule that exponents would only auto-combine if
they had the same term minus the coefficient, the same rule that `Add` uses
to determine what terms should auto combined by addition. So it would
combine `exp(x)*exp(x)` into `e**(2x)$, but `exp(x)*exp(y)` would be left
alone. It turns out that some of our algorithms, namely the Gruntz limit
algorithm, relies on auto-combining. We already had a function that could
combine exponents, `powsimp()`, but it also combined bases, as in
`x**z*y**z` into `(x*y)**z`, so I had to split the behavior so that it
could act only as auto-combining once did (by the way, use `powsimp(expr,
combine='exp', deep=True)` to do this). Then, after some help from Ondrej
on pinpointing the exact location of the bugs, I just applied the function
there. The last thing I did here was to split the behavior of expand, so
that you could do`expand(x*(y + 1), mul=False)` and it would leave it
alone, but`expand(exp(x + y), mul=False)` would return `exp(x)*exp(y)`.
This split behavior turned out to be useful in more than one place in my
code.
-
-This was the first non bug fix patch of mine that was pushed in to SymPy,
and at the time of this writing, it is the last major one in the latest
stable version. It took some major rebasing to get my convoluted commit
history ready for submission, and it was during this phase that I git
finally clicked for me, especial the <code>git rebase</code> command. This
work took a few weeks from my ODEs time, and it became clear that I would
not be doing every possible thing from my application. The reason that I
included so much in my application was that my project was non-atomic. I
could implement a little or a lot and still have a working useful module.
+So I spent a few weeks hacking on the core to make it not auto-combine
exponents. I came up with a rule that exponents would only auto-combine if
they had the same term minus the coefficient, the same rule that `Add` uses
to determine what terms should auto combined by addition. So it would
combine `exp(x)*exp(x)` into `e**(2x)`, but `exp(x)*exp(y)` would be left
alone. It turns out that some of our algorithms, namely the Gruntz limit
algorithm, relies on auto-combining. We already had a function that could
combine exponents, `powsimp()`, but it also combined bases, as in
`x**z*y**z` into `(x*y)**z`, so I had to split the behavior so that it
could act only as auto-combining once did (by the way, use `powsimp(expr,
combine='exp', deep=True)` to do this). Then, after some help from Ondrej
on pinpointing the exact location of the bugs, I just applied the function
there. The last thing I did here was to split the behavior of expand, so
that you could do`expand(x*(y + 1), mul=False)` and it would leave it
alone, but`expand(exp(x + y), mul=False)` would return `exp(x)*exp(y)`.
This split behavior turned out to be useful in more than one place in my
code.
+
+This was the first non bug fix patch of mine that was pushed in to SymPy,
and at the time of this writing, it is the last major one in the latest
stable version. It took some major rebasing to get my convoluted commit
history ready for submission, and it was during this phase that I git
finally clicked for me, especially the `git rebase` command. This work
took a few weeks from my ODEs time, and it became clear that I would not be
doing every possible thing from my application. The reason that I included
so much in my application was that my project was non-atomic. I could
implement a little or a lot and still have a working useful module.
If you look at my timeline on my application, you can see that the first
half is symbolic methods, and the second half is other methods, things like
series. It turns out that we didn't really learn much about systems of
ODEs in my course and we learned very little about numerical methods (and
it would take much more to know how to implement them). We did learn
series methods, but they were so annoying to do that I came to hate them
with a passion. So I decided to just focus on symbolic methods, which were
my favorite anyway. My goal was to implement as many as I could.
-After I finished up separable equations, I came up with an idea that I did
not have during the application period. <code>dsolve()</code> was becoming
cluttered fast with all of my solution methods. The way that it worked was
that it took an ODE and it tried to match methods one by one until it found
one that worked, which it then used. This had some drawbacks. First, as
I mentioned, the code was very cluttered. Second, the ODEs methods would
have to be applied in a predetermined order. There are several ODEs that
match more than one method. For example, `2*x*y + (x**2 + y**2)*dy/dx=0$
has coefficients that are both homogeneous of order 2, and is also exact,
so it can be solved by either method. The two solvers return differently
formatted solutions for each one. A simpler example is that 1st order ODEs
with homogeneous coefficients can be solved in two different ways. My
working solution was to try them both and then apply some heuristics to
return the simplest one. But sometimes, one way would create an impossible
integral that would hand the integration engine. And it made debugging the
two solvers more difficult because I had to override my heuristic. This
also touches on the third point. Sometimes the solution to an ODE can only
be represented in the form of an unevaluatable integral. SymPy's
<code>integrate()</code> function is supposed to return an unevaluated
<code>Integral</code> class if it cannot do it, but all too often it will
just hang forever.
-
-The solution I came up with was to rewrite dsolve using a hints method. I
would create a new function called <code>classify_ode()</code> that would
do all of the ODE classification, removing it from the solving code. By
default, dsolve would still use a predetermined order of matching methods.
But you could override it by passing a "hint" to <code>dsolve</code> for
any matching method, and it would apply that method. There would also be
options to only return unevaluated integrals when applicable.
-
-I ended up doing this and more (see the docstrings for
<code>classify_ode()</code> and <code>dsolve()</code> in the current git
master branch), but before I could I needed to clean up some things. I
needed to rewrite all of <code>dsolve()</code> and related functions.
Before I started the program, there were some special cases in dsolve for
second order linear homogeneous ODEs with constant coefficients and one
very special case ODE for the expanded form of $latex
\frac{d^2}{dx^2}(xe^{-y}) = 0$.
-
-So the first thing I did was implement a solver for the general
homogeneous linear with constant coefficients case. These are rather
simple to do: you just find the roots of the characteristic polynomial
built off of the coefficients, and then put the real parts of the roots in
front of the argument of an exponential and the imaginary parts in front of
the arguments of a sine and cosine (for example, $latex 3 \pm 2i$ would
give $latex C1e^{3x}\sin{2x} + C2e^{3x}\cos{2x}$. The thing was, that if
the imaginary part is 0, then you only have 1 arbitrary constant on the
exponential, but if it is non-zero, you get 2, one for each trig function.
The rest falls out nicely if you plug 0 in for $latex b$ into
$e^{ax}(C1\sin{bx} + C2\cos{box})$ because the sine goes to 0 and the
cosine becomes 1. But you would end up with $latex C1 + C2$ instead of
just $latex C1$ in that case. I had already planned on doing arbitrary
constant simplification as part of my project, so I figured I would put
this on hold and do that first. Then, once that was done, the homogeneous
case would be reduced to 1 case instead of the usual 2 or 3.
-
-My original plan was to make an arbitrary constant type that automatically
simplified itself. So, for example, if you entered <code>C1 + 2 + x</code>
with <code>C1</code> an arbitrary constant, it would reduce to just
<code>C1 + x</code>. I worked with Ondrej, including visiting him in Los
Alamos, and we build up a class that worked. The problem was that, in
order to have auto-simplification, I had to write the simplification
directly into the core. Neither of us liked this, so we worked a little
bit on a basic core that would allow auto-simplification to be written
directly in the classes instead of in the <code>Mul.flatten()</code> and
<code>Add.flatten()</code> methods. It turns out that my constant class
isn't the only thing that would benefit from this. Things like the order
class (O(x)) and the infinity class (oo) are auto-simplified in the core,
and things could be much cleaner if they happened in the classes
themselves. Unfortunately, modifying the core like this is not something
that can happen overnight or even in a few weeks. For one thing, it needed
to wait until we had the new assumptions system, which was another Google
Summer of Code project running parallel to my own. So we decided to shelf
the idea.
-
-I still wanted constant simplification, so I settled with writing a
function that could do it instead. There were some downsides to this.
Making the function as general as the classes might have been would have
been far too much work, so I settled on making it an internal-only
function that only worked on symbols named <code>C1</code>,
<code>C2</code>, etc. Also, unlike writing the simplification straight
into <code>Mul.flatten()</code> which was as simple as removing any terms
that were not dependent on x, writing a function that parsed an expression
and simplified it was considerably harder to write. I managed to churn out
something that worked, and so I was ready to finish up the solver I had
started a few paragraphs ago.
-
-After I finished that, I still needed to maintain the ability to solve
that special case ODE. Apparently, it is an ODE that you would get
somewhere in deriving something about relativity, because it was in the
relativity.py example file. I used Maple's excellent
<code>odeanalyser()</code> function (this is where I go the idea for my
<code>classify_ode()</code>)to find a simple general case ODE that it fit
(Liouville ODE). After I finished this, I was ready to start working on
the hints engine.
-
-It took me about a week to move all classification code into
<code>classify_ode()</code>, move all solvers into individual functions,
separate simplification code into yet other functions, and tie it all
together in <code>dsolve()</code>. In the end, the model worked very
well. The modularization allowed me to do some other things that I had not
considered, such as creating a special "best" hint that used some
heuristics that I originally developed for first order homogeneous which
always has two possible solutions to try to give the best formatted
solution for any ODE that has more than one possible solution method. It
also made debugging individual methods much easier, because I could just
use the built in hint calls in <code>dsolve()</code> instead of commenting
out lines of code in the source.
+After I finished up separable equations, I came up with an idea that I did
not have during the application period. `dsolve()` was becoming cluttered
fast with all of my solution methods. The way that it worked was that it
took an ODE and it tried to match methods one by one until it found one
that worked, which it then used. This had some drawbacks. First, as I
mentioned, the code was very cluttered. Second, the ODEs methods would
have to be applied in a predetermined order. There are several ODEs that
match more than one method. For example, `2*x*y + (x**2 + y**2)*dy/dx=0`
has coefficients that are both homogeneous of order 2, and is also exact,
so it can be solved by either method. The two solvers return differently
formatted solutions for each one. A simpler example is that 1st order ODEs
with homogeneous coefficients can be solved in two different ways. My
working solution was to try them both and then apply some heuristics to
return the simplest one. But sometimes, one way would create an impossible
integral that would hand the integration engine. And it made debugging the
two solvers more difficult because I had to override my heuristic. This
also touches on the third point. Sometimes the solution to an ODE can only
be represented in the form of an unevaluatable integral. SymPy's
`integrate()` function is supposed to return an unevaluated `Integral`
class if it cannot do it, but all too often it will just hang forever.
+
+The solution I came up with was to rewrite dsolve using a hints method. I
would create a new function called `classify_ode()` that would do all of
the ODE classification, removing it from the solving code. By default,
dsolve would still use a predetermined order of matching methods. But you
could override it by passing a "hint" to `dsolve` for any matching method,
and it would apply that method. There would also be options to only return
unevaluated integrals when applicable.
+
+I ended up doing this and more (see the docstrings for `classify_ode()`
and `dsolve()` in the current git master branch), but before I could I
needed to clean up some things. I needed to rewrite all of `dsolve()` and
related functions. Before I started the program, there were some special
cases in dsolve for second order linear homogeneous ODEs with constant
coefficients and one very special case ODE for the expanded form of
`d^2/dx^2 x*e**(-y) = 0`.
+
+So the first thing I did was implement a solver for the general
homogeneous linear with constant coefficients case. These are rather
simple to do: you just find the roots of the characteristic polynomial
built off of the coefficients, and then put the real parts of the roots in
front of the argument of an exponential and the imaginary parts in front of
the arguments of a sine and cosine (for example, `3 +/- 2i` would give
`C1*e**(3*x)*sin(2*x) + C2*e**(3*x)*cos(2*x)`. The thing was, that if the
imaginary part is 0, then you only have 1 arbitrary constant on the
exponential, but if it is non-zero, you get 2, one for each trig function.
The rest falls out nicely if you plug 0 in for `b` into
`e**(ax)*(C1*sin(b*x) + C2*cos(b*x)` because the sine goes to 0 and the
cosine becomes 1. But you would end up with`C1 + C2` instead of just `C1`
in that case. I had already planned on doing arbitrary constant
simplification as part of my project, so I figured I would put this on hold
and do that first. Then, once that was done, the homogeneous case would be
reduced to 1 case instead of the usual 2 or 3.
+
+My original plan was to make an arbitrary constant type that automatically
simplified itself. So, for example, if you entered `C1 + 2 + x` with `C1`
an arbitrary constant, it would reduce to just `C1 + x`. I worked with
Ondrej, including visiting him in Los Alamos, and we build up a class that
worked. The problem was that, in order to have auto-simplification, I had
to write the simplification directly into the core. Neither of us liked
this, so we worked a little bit on a basic core that would allow
auto-simplification to be written directly in the classes instead of in the
`Mul.flatten()`and `Add.flatten()` methods. It turns out that my constant
class isn't the only thing that would benefit from this. Things like the
order class (`O(x)`) and the infinity class (`oo`) are auto-simplified in
the core, and things could be much cleaner if they happened in the classes
themselves. Unfortunately, modifying the core like this is not something
that can happen overnight or even in a few weeks. For one thing, it needed
to wait until we had the new assumptions system, which was another Google
Summer of Code project running parallel to my own. So we decided to shelf
the idea.
+
+I still wanted constant simplification, so I settled with writing a
function that could do it instead. There were some downsides to this.
Making the function as general as the classes might have been would have
been far too much work, so I settled on making it an internal-only
function that only worked on symbols named `C1`, `C2`, etc. Also, unlike
writing the simplification straight into `Mul.flatten()` which was as
simple as removing any terms that were not dependent on x, writing a
function that parsed an expression and simplified it was considerably
harder to write. I managed to churn out something that worked, and so I
was ready to finish up the solver I had started a few paragraphs ago.
+
+After I finished that, I still needed to maintain the ability to solve
that special case ODE. Apparently, it is an ODE that you would get
somewhere in deriving something about relativity, because it was in the
relativity.py example file. I used Maple's excellent `odeanalyser()`
function (this is where I go the idea for my `classify_ode()`)to find a
simple general case ODE that it fit (Liouville ODE). After I finished
this, I was ready to start working on the hints engine.
+
+It took me about a week to move all classification code into
`classify_ode()`, move all solvers into individual functions, separate
simplification code into yet other functions, and tie it all together in
`dsolve()`. In the end, the model worked very well. The modularization
allowed me to do some other things that I had not considered, such as
creating a special "best" hint that used some heuristics that I originally
developed for first order homogeneous which always has two possible
solutions to try to give the best formatted solution for any ODE that has
more than one possible solution method. It also made debugging individual
methods much easier, because I could just use the built in hint calls in
`dsolve()` instead of commenting out lines of code in the source.
This was good, because there was one more method that I wanted to
implement. I wanted to be able to solve the inhomogeneous case of a nth
order linear ode with constant coefficients. This can be done in the
general case using the method of variation of parameters. It was quite
simple to set up variation of parameters up in the code. You only have to
set up a system of integrals using the Wronskian of the general solutions.
It would usually be a very poor choice of a method if you were trying to
solve an ODE by hand because taking the Wronskian and computing n integrals
is a lot of work. But for a CAS, the work is already there. I just have
to set up the integrals.
@@ -89,3 +89,6 @@
== Post-GSoC ==
I plan on continuing development with SymPy now that the Google Summer of
Code period is over. SymPy is an amazing project, mixing Python and
Computer Algebra, and I want to help it grow. I may even apply again in a
future year to implement some other thing in SymPy, or maybe apply as a
mentor for SymPy to help someone else improve it.
+
+== Commits ==
+All of my commits on
[http://github.com/asmeurer/sympy/commits/odes-GSoC-2009 this branch] fell
under the 2009 Google Summer of Code period
(19b33d943e06a3b4ac9cdebb41901af853dc2224).