On Mon, Aug 25, 2008 at 12:55 AM, Jason Merrill <[EMAIL PROTECTED]> wrote:
>
> In hopes that it may be a useful reference during the current work on
> symbolics, I wrote a toy Mathematica program for transforming a single
> higher order ODE into a system of first order ODEs. Most of the free
> numerical differential equation solvers I've seen want input in the
> form y'[x] == f(y,x), so the purpose of a program like this is to take
> more general input and turn it into something suitable for those
> routines.
>
> I'll use the Poisson-Boltzmann equation in spherical coordinates, a
> second order ODE, as my test example, because I recently needed to
> solve it numerically.
>
> pbe = r^2 p''[r] + 2 r p'[r] == r^2 k^2 Sinh[p[r]]
>
> The full form of the equation, which is useful for pattern matching,
> is
>
> FullForm[pbe]
>
> Equal[Plus[Times[2,r,Derivative[1][p][r]],Times[Power[r,
> 2],Derivative[2][p][r]]],Times[Power[k,2],Power[r,2],Sinh[p[r]]]]
>
> Notice how this looks a lot like Lisp?
>
> First, I want to find out what the order of the equation is, so I will
> match any terms of the form Derivative[o][p][r], and then find the
> match with the maximum o:
>
> diffEqOrder[eqn_, y_, x_] :=
> Max[Cases[pbe, Derivative[o_][y][x] -> o, Infinity]]
>
> Cases produces a list of all the subexpressions in its first argument
> that match the pattern specified in the second argument, optionally
> with a replacement rule applied to the matches. By default, it only
> looks at the top level of an expression, but setting the final
> argument to Infinity forces it to look all the way down the tree. The
> _ here is a wildcard that matches any expression, and writing o_ gives
> that subexpression a name. Since o is the only part I care about, I
> replace all the matches with just o.
>
> diffEqOrder[pbe, p, r]
> 2
>
> Now we're ready for the function that actually transforms equations:
>
> firstOrderForm[eqn_, y_, x_] :=
> Module[{rep, order = diffEqOrder[eqn, y, x]},
> rep = Solve[
> eqn /. {Derivative[o_][y][x] -> y[o][x], y[x] -> y[0][x]},
> y[order][x]];
> Join[Table[y'[o][x] == y[o + 1][x], {o, 0, order - 2}],
> y'[order - 1][x] == y[order][x] /. rep]
> ]
>
> The Module bit is just for keeping variables used in the definition
> out of global scope. The important line is
>
> eqn /. {Derivative[o_][y][x] -> y[o][x], y[x] -> y[0][x]}
>
> where I use the replacement operator /. to replace all the derivatives
> of y with respect to x by a new function of x, y[o][x], where o is
> again the order of the derivative. Once I've done that, I solve the
> equation for the highest order term. Then I make a list basically
> just giving definitions of all the new functions as derivatives of the
> next lower order functions, and stick the solved highest order
> derivative equation on the end of the list.
>
> Here's the result for the example equation
>
> firstOrderForm[pbe, p, r]
>
> {p'[0][r] == p[1][r], p'[1][r] == (k^2 r Sinh[p[0][r]] - 2 p[1][r])/r}
>
> It's just a short shot from here to providing the input in a form that
> would be accepted by the gsl solver, or scipy.
>
> Hopefully this will be useful as one example of the kind of
> manipulations that pattern matching makes possible. I'm happy to
> answer any more questions about it.
Thanks for doing this. It is indeed very useful, I always wondered how
things like this are done in Mathematica.
Any ideas how this could be nicely translated to Python?
Ondrej
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---
