Separable next

May 31, 2009

So the next thing on my list is seperable differential equations. As I predicted in my proposal, I am going to have to do some work to get SymPy to recognize all equations that are separable. Currently, it can match expressions that are exactly written as $a(x)b(y)\frac{dy}{dx}+c(x)d(y)=0$, such as $x^2(y+2)\frac{dy}{dx}-y^3=0$ but the matching engine cannot get things like $(x^2y+2x^2)\frac{dy}{dx}+-y^3=0$, because it doesn’t know how to factor out the $x^2$. There is a function, collect, that can factor out expressions if you give them to it explicily with

collect(x**2*y+x**2,x**2)

but it doesn’t know yet how to separate variables in general. So I will be working on it this week. Once it can properly separate separable expressions, implementing it in dsolve will be cake.

Also, it needs to learn how to do $e^{x+y}\rightarrow e^{x}e^{y}$ and $1+x+y+xy \rightarrow (1+x)(1+y)$.

First Order Differential Equations with Homogeneous Coefficients

May 31, 2009

So here is what I have done.  Homogeneous has several meanings, but when I use it in this post, I will mean homogeneous functions.  We only care about functions of two variables, so I will also only talk about them.  A function $F(x,y)$ is called homogeneous of order $n$ if $F(tx,ty)=t^nF(x,y)$ for some $n$.  Also, if we “pull” an $x$ out, we get $F(x,y)=x^nF(1,\frac{x}{y})=x^nG(\frac{y}{x})$, which is why the definition is often stated that the function can be written as a function of $\frac{y}{x}$.  Of course, we could do the same with $y$ and get a function of $\frac{x}{y}$.  Here is an example of a homogeneous function: $x^2+y\sqrt{x^2+y^2}\sin{\frac{x}{y}}$.  You can see that $(xt)^2+yt\sqrt{(xt)^2+(yt)^2}\sin{\frac{xt}{yt}}=t^2\left(x^2+y\sqrt{x^2+y^2}\sin{\frac{x}{y}}\right)$.

Now, if the coefficients $P$ and $Q$ in the first order differential equation $P(x,y)dx+Q(x,y)dy=0$ are both homogeneous functions of the same order, then the substitution $x=uy$ or $y=ux$ will make the equation separable.  I will show how to do it for $y=ux$.  We have $dy=xdy+udx$ by the product rule.  We also have $P(x,y)=P(x,xu)$ and $Q(x,y)=Q(x,xu)$.  Because these functions are homogeneous, we can “pull out” an $x^n$ from each, giving us $x^nP(1,u)dx+x^nQ(1,u)(xdu+udx)=0$.  Because we made the substitution $u=\frac{y}{x}$, we already have had to assume $x\neq0$.  Thus, we can divide the whole thing by $x^n$.  Doing this, and expanding the remaining terms, we get $P(1,u)dx+xQ(1,u)du+uQ(1,u)dx=0$.  Refactoring the $dx$ and $dy$ terms, we get $(P(1,u)+uQ(1,u))dx=-xQ(1,u)du$.  The equation is separable!  Separating variables and integrating, we get $\int{\frac{dx}{x}}=\int{\frac{-Q(1,u)du}{P(1,u)+uQ(1,u)}}+C$.  If we had used $x=uy$ instead, we would have gotten $\int{\frac{dy}{y}}=\int{\frac{-P(u,1)du}{uP(u,1)+Q(u,1)}}+C$.  This is exactly how I was able to solve these equations in SymPy.  Each homogeneous equation has two possible integrals, and often the right hand side of one equation is a much harder integral than the right hand side of the other.  Therefore, I did them both, and applied a little heuristic on which one to return.  It prefers expressions that can be solved for y (f(x) in the case of SymPy), expressions that have evaluated integrals, and if neither of them are solvable, the shortest one, which tends to be the simplest.

Here is an example from my text book.  $2xydx+(x^2+y^2)dy=0$.  Astute readers of this blog may have noticed that this equation is exact.  It is easier to solve that way, but we will try using the methods outlined above.  First, we notice that the coeficients of $dx$ and $dy$ are both homogeneous of order 2, so we can make the substitution $x=uy$ or $u=yx$ and the equation would become seperable.  If I were doing this by hand or for a homework assignment, I would make one of those substitutions and work it through, but we have the exact form above, so let’s use it just like SymPy would.  $P(x,y)=2xy$ and $Q(x,y)=x^2+y^2$.  So we get $\int{\frac{dx}{x}}=\int{\frac{-(1+u^2)du}{2u+u(1+u^2)}}+C$ with $u=\frac{y}{x}$ and $\int{\frac{dy}{y}}=\int{\frac{-2udu}{u(2u)+(1+u^2)}}+C$ with $u=\frac{x}{y}$ (here and for the remainder of this derivation, the arbitrary constants in the two equations do not necessarily equal each other until the last step).  Both integrals can be solved  with a substitution of the denominator, giving, with $u$ back substituted, $\ln{x}=\frac{-\ln{(\frac{y^3}{x^3}+\frac{3y}{x})}}{3}+C$ and $\ln{y}=\frac{-\ln{(\frac{3x^3}{y^3}+1)}}{3}+C$.  If we make $C=ln(C_1)$, and pull the constants in front of the logs on the right hand sides in as powers, we can combine everything into logarithms.  After doing that, we then take the antilog of both sides and get $x=C_1{(\frac{y^3}{x^3}+\frac{3y}{x})}^{-\frac{1}{3}}$ and $y=C_1{(\frac{3x^3}{y^3}+1)}^{-\frac{1}{3}}$.  We then raise the whole equation to the $-3$ power and make $C_1^{-3}=A$.  We get
$\frac{1}{x^3}=A(\frac{y^3}{x^3}+3\frac{y}{x})$
$\frac{1}{x^3}=A\frac{y^3+3x^2y}{x^3}$
$K=y^3+3x^2y$
and
$\frac{1}{y^3}=A(\frac{3x^2}{y^2}+1)$
$\frac{1}{y}=A(3x^2y+y^2)$
$K=y^3+3x^2y$,
where $K=\frac{1}{A}$

Verify as an exercise that you get the same solution by the exact differential equation methods (or don’t… you know, depending on which you’d rather do).

Of course, if that was all there was to it, I would have finished in a day.  There was more.  First, I had to write a function that determines if an expression is homogeneous and what order if it is. That was about 150 lines of code right there plus tests (150 lines is a lot in Python).

I also wanted it to recognize that $\ln{x}-\ln{y}$ is homogeneous because $\ln{x}-\ln{y}=\ln{\frac{x}{y}}$. SymPy was incapable of combining logarithms like this, so I had to write a logcombine function, which was another 100 lines of code plus tests. This came in handy, because I was also able to use it in the part above where I combined the logarithms of my answer to eliminate all of them.

Lastly, you may have noticed that my math looks much nicer in this post. It turns out that WordPress supports math using \$ Read the rest of this entry »

Just finished a new bit.

May 30, 2009

I just finished first order homogeneous equations, as well as my logcombine function, and I have pushed them to http://github.com/asmeurer/sympy/tree/odes-stable.  So expect a major blog posting about them both soon (it is midnight right now, so bedtime, maybe tomorrow).

If you are a member of the SymPy project reading this, please, review my patches.  I have already submitted exact equations and the logcombine function, and neither of them have made it in.  And the more reviewing that there is, the better for my code.

One problem solved

May 25, 2009

Well, I figured out why my patch files are appearing corrupted to everyone.  It is close to what I expected.  The files had a Macintosh character encoding, when everyone was expecting UTF-8.  Changing the encoding with TextWrangler was easy enough, but now the question is, how do I make git use this encoding by default?

May 25, 2009

So the GSoC program officially started Saturday.  I have almost finished everything on the timeline for the first week already, except for initial conditions, which I am putting off until Ondrej recodes solve so that it returns a RootOf instance when it cannot solve the equation.  See this post.

In order to solve homogeneous equations, I had to determine if an equation is homogeneous.  I will post what all this means, when I finish everything up.  Right now, I just need to write a few more tests and make sure that they work.  So I wrote the function, which is very complicated, because it has to take any expression and break it up recursively.  One thing that I wanted for it to do is to realize that log(x)-log(y) == log(x/y) (no LaTeX today, I am too tired).  This is basically because just log(x) or log(y) will make it non-homogeneous, but log(x/y) will make is homogeneous (again, more on what this all means later).  It turns out SymPy could not really do this simplification, so I had to go in and write a logcombine function that did that.  The function works quite well (at least from my tests), but it is also a beefy piece or recursive work.

So here is the part where I am now.  I wanted to submit the logcombine function as a regular patch, because strictly speaking it is not related to my GSoC project.  So I used ‘git format-patch’ and git ‘send-email’ as usual.  But then, yesterday, I began to wonder why no one had reviewed or even responded to my patches.  I checked that patch list, and found to my dismay that none of the emails made it through.  I suspect that git is the problem, so I send the patches manually with my email client.

Which led to another problem: the reviewers are claiming that my files are corrupt.  These are major problems.  I cannot reliably send patch emails through git, and when I send them manually, the files are messed up.

Can anyone reading this blog (does anyone actually read this blog?) give me some advice here?  I think the corruption of the patch files may have something to do with the fact that I am using Mac OS X, which I know likes to put nasty things like resource forks in files which screws them up when they are sent to non-Macs.

Of course, this does not explain why git will not send emails.  It has sent them before, though I think I have had this problem since the beginning.

And another question to you git geniuses out there: how do I change branches without committing changes?  Sometimes, I am in the middle of two things, and I finish one and I want to push it to its own branch, but I cannot do this without committing the changes on the other thing, which I am not ready to do because I am not done with it yet.

Work started — Exact Differential Equations

May 16, 2009

I have decided to start working on my project, even though it is a week early, since I am finished with school and have nothing better to do.  I also like the idea of getting ahead.  I am pushing everything to my github account.  So far, I have implemented exact differential equations, which are cake (i.e., they took me about an hour to code up).  These are equations of the type,where.  If this condition holds, then there exists a function, often called a ‘potential function’ because of some applications of theses equations to physics, such thatand.  This is because mixed partials are equal for continuous functions.  The solution will then just be.  The tricky part is finding the potential function, but it turns out to be easier than you would think.  Because of the fundamental theorem of calculus, the potential function is just.  There is a restriction where the rectangle connecting,,, andhas to lie entirely in the domain of and, but if we let, then this is not really a problem for functions that SymPy will encounter.  UPDATE: It turns out this is a problem if the equation has singularities in it.  Fortunately, I was able to code up a workaround that usually works.

So you can see that this is easy to implement in SymPy if you know the above fact.  It turns out that most solution methods for ODEs are like this.  They can be solved in the general case, although usually students are only taught tricks because they are much easier to remember than generally solved formulas.  For example, to solve an exact differential equation, students are often taught to just integratewith respect toandwith respect to.  It is much simpler for a human being to do that than the above integral, because the integral involves evaluating limits and so on, but for a computer algebra system, the above integral is a one-liner.

By the way, if you want to test my code, you should clone my github repository and switch to the odes branch.  You can do this by typing ‘git clone git://github.com/asmeurer/sympy.git’ and then ‘git checkout odes’ in a Terminal that is cd’d to the directory you want to clone to (of course, you will need git installed first!).  Then type ‘cd sympy’ and ‘./bin/isympy’, which will start SymPy.  Then from there, you can do stuff like ‘dsolve(sin(x)*cos(f(x))+cos(x)*sin(f(x))*f(x).diff(x),f(x)’ (this is exact).  It’s easy to generate an exact differential equation.  Just start with a random function of x and y, then take the derivative of it with respect to x and y each, and do a subs(y,f(x)) on each term (SymPy wants functions for the dependent variable, so we have to use f(x)).  Then do dsolve(<the derivative with respect to x>+<the derivative with respect to y>*f(x).diff(x),f(x)).  You should get your original function equals an arbitrary constant C1.  By the way, diff(<function>,x) will take the derivative of <function> with respect to x.

Up next: Initial conditions, and then first order homogeneous differential equations.

My GSoC Proposal

May 8, 2009

I have copied my proposals over to the SymPy Wiki and formatted them so that they look nice (including LaTeX for the formulas).  I included both the Portland State University proposal, which was accepted (at the top), and the Python Software Foundation proposal, which was nearly identical, but was not accepted.  I did this mainly for myself so I can easily find my timeline during the summer, but if anyone wants to see what my proposal looked like, here it is.  Note that I am pretty sure that the timeline on the proposal will not be the same as what I will actually do (e.g., I have learned much less about systems of ODEs (0), than I expected, and I have learned about a couple new symbolic methods).