On Aug 24, 6:55 pm, 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.
>
> Regards,
>
> Jason Merrill
I'm afraid I'm straying off topic, but I thought I'd use the above
code as a jumping off point to say a few more things about pattern
matching and Mathematica, since there seems to be at least some
interesting in hearing about this kind of thing. First, I should
point out that there was a bug in my definition of diffEqOrder:
> diffEqOrder[eqn_, y_, x_] :=
> Max[Cases[pbe, Derivative[o_][y][x] -> o, Infinity]]
should be
diffEqOrder[eqn_, y_, x_] :=
Max[Cases[eqn, Derivative[o_][y][x] -> o, Infinity]] (*Changed pbe
to eqn on this line*)
Pattern matching is one of the fundamental underpinnings of
Mathematica--more than is obvious at first pass. I used it
conspicuously to play with derivative terms, but it is used
inconspicuously in every function definition. Remember how _ is a
wildcard, and o_ gives the match a name so it can be used again
later? Well that's what's going on in function definitions. In its
definition, each of the arguments of diffEqOrder is a named pattern,
which is why we can use the names eqn, y, and x later in the
definition. So function definitions are just replacement rules that
are applied globally after they are defined.
You can be more restrictive with the argument patterns than just using
wildcards. For instance, say I wanted to allow the first argument to
just be an expression, with an implied == 0 at the end. I could have
done
Clear[diffEqOrder] (*just gets rid of the old definition*)
diffEqOrder[eqn_Equal,y_,x_] :=
Max[Cases[eqn, Derivative[o_][y][x] -> o, Infinity]]
diffEqOrder[expr_,y_,x_] := diffEqOrder[expr == 0, y, x]
The eqn_Equal is a more restrictive pattern that only matches
expressions whose head is Equal, i.e. Equal[_,_] which is the same as
_ == _ . Here's a factorial definition in the same style
fac[n_Integer] := n*fac[n - 1]
fac[1] := 1
Part of what allows pattern matching to be so powerful is that
*everything* has a full form that is just an expression tree tree,
like
> 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]]]]
Even plots have this form
FullForm[Plot[Sin[x], {x, -Pi, Pi}]]
Graphics[List[List[List[],List[],List[Hue[0.67`,0.6`,
0.6`],Line[List[List[-3.1415925253615216`,-1.2822827167891634`*^-7],
<...>,List[3.1415925253615216`,
1.2822827167891634`*^-7]]]]]],List[Rule[AspectRatio,Power[GoldenRatio,-1]],Rule[Axes,True],Rule[AxesOrigin,List[0,0]],
Rule[PlotRange,List[List[Times[-1,Pi],Pi],List[-0.9999998782112116`,
0.9999998592131705`]]],
Rule[PlotRangeClipping,True],Rule[PlotRangePadding,List[Scaled[0.02`],Scaled[0.02`]]]]]
Not that friendly to look at, but this form lets me manipulate plots
just like any other data. If I want to see just the points that are
evaluated for the plot, omitting the line connecting them, I can do
Plot[Sin[x],{x,-Pi,Pi}] /. Line->Point
This gets to the reason that I think Mathematica is actually a
beautiful language. First, it has the Lisp thing going for it--
everything is an expression, and programs are data. Given that
foundation, it uses pattern matching to do basically all of the work
with expressions, which makes things feel really unified once you get
it. It also provides powerful functional programming constructs: Map,
Nest, Fold, Thread, and variations on these themes.
Mathematica is nominally "multi-paradigm": it has For loops and other
procedural stuff; it also has objects, sort of, though I've never seen
much done with them. But really, functional programming and pattern
matching are the beautiful part, and they let you do everything. If
you've tried to program Mathematica a different way, you probably
thought it was ugly. But then, you were doing it wrong.
Cheers,
JM
--~--~---------~--~----~------------~-------~--~----~
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
-~----------~----~----~----~------~----~------~--~---
