* algebra/odeef.spad.pamphlet (ElementaryFunctionODESolver): Tidy.

1 files changed, 26 insertions, 11 deletions

diff --git a/src/algebra/odeef.spad.pamphlet b/src/algebra/odeef.spad.pamphlet
index 5598124b..37f26e8e 100644
--- a/src/algebra/odeef.spad.pamphlet
+++ b/src/algebra/odeef.spad.pamphlet
@@ -434,7 +434,7 @@ ElementaryFunctionODESolver(R, F): Exports == Implementation where
     getfreelincoeff1: (F, K, List F) -> F
     getlincoeff     : (F, K) -> F
     getcoeff        : (F, K) -> UU
-    parseODE        : (F, OP, SY) -> Union(LEQ, NLQ)
+    parseODE        : (F, OP, SY) -> Union(LEQ, NLQ, "failed")
     parseLODE       : (F, List K, UP, SY) -> LEQ
     parseSYS        : (List F, List OP, SY) -> Union(SYS, "failed")
     parseSYSeq      : (F, List K, List K, List F, SY) -> Union(ROW, "failed")
@@ -462,7 +462,9 @@ ElementaryFunctionODESolver(R, F): Exports == Implementation where
     solve(diffeq:F, y:OP, center:EQ, y0:List F) ==
       a := rhs center
       kx:K := kernel(x := retract(lhs(center))@SY)
-      (ur := parseODE(diffeq, y, x)) case NLQ =>
+      ur := parseODE(diffeq, y, x)
+      ur case "failed" => "failed"
+      ur case NLQ =>
         not one?(#y0) => error "solve: more than one initial condition!"
         rc := ur::NLQ
         (u := solve(rc.dx, rc.dy, y, x)) case "failed" => "failed"
@@ -483,7 +485,9 @@ ElementaryFunctionODESolver(R, F): Exports == Implementation where
         error "solve: not a first order linear system"
 
     solve(diffeq:F, y:OP, x:SY) ==
-      (u := parseODE(diffeq, y, x)) case NLQ =>
+      u := parseODE(diffeq, y, x)
+      u case "failed" => "failed"
+      u case NLQ =>
         rc := u::NLQ
         (uu := solve(rc.dx, rc.dy, y, x)) case "failed" => "failed"
         uu::F
@@ -527,28 +531,39 @@ ElementaryFunctionODESolver(R, F): Exports == Implementation where
           eq := eq - ci * y::F
       [n, v, -eq]
 
--- returns either [p, g] where the equation (diffeq) is of the form p(D)(y) = g
--- or [p, q] such that the equation (diffeq) is of the form p dx + q dy = 0
+    -- returns either 
+    --   [p, g] where the equation (diffeq) is of the form 
+    --        p(D)(y) = g
+    -- or
+    --   [p, q] such that the equation (diffeq) is of the form 
+    --        p dx + q dy = 0
     parseODE(diffeq, y, x) ==
-      f := y(x::F)
+      -- FIXME: Assume we don't have a parametric equation.
+      f := y(x::F)   
       l:List(K) := [retract(f)@K]
       n:N := 2
       for k in varselect(kernels diffeq, x) | is?(k, OPDIFF) repeat
         if (m := height k) > n then n := m
       n := (n - 2)::N
--- build a list of kernels in the order [y^(n)(x),...,y''(x),y'(x),y(x)]
+      -- build a list of kernels in the order [y^(n)(x),...,y''(x),y'(x),y(x)]
       for i in 1..n repeat
         l := concat(retract(f := differentiate(f, x))@K, l)
+      -- Determine the order of the equation and its leading coefficient
       k:K   -- #$^#& compiler requires this line and the next one too...
       c:F
-      while not(empty? l) and zero?(c := getlincoeff(diffeq, k := first l))
-        repeat l := rest l
-      empty? l or empty? rest l => error "parseODE: equation has order 0"
+      while not empty? l repeat 
+        c' := getcoeff(diffeq, k := first l)
+        c' case "failed" => return "failed"
+        c := c'::F
+        not zero? c => break
+        l := rest l
+      empty? l or empty? rest l => "failed" -- maybe a functional equation?
       diffeq := diffeq - c * (k::F)
       ny := name y
       l := rest l
       height(k) > 3 => parseLODE(diffeq, l, monomial(c, #l), ny)
-      (u := getcoeff(diffeq, k := first l)) case "failed" => [diffeq, c]
+      u := getcoeff(diffeq, k := first l)
+      u case "failed" => [diffeq, c]
       eqrhs := (d := u::F) * (k::F) - diffeq
       freeOf?(eqrhs, ny) and freeOf?(c, ny) and freeOf?(d, ny) =>
         [monomial(c, 1) + d::UP, eqrhs]