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.
Regards,
Jason Merrill
--~--~---------~--~----~------------~-------~--~----~
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
-~----------~----~----~----~------~----~------~--~---
