I have been extending the ode2() function to solve some non-linear first
order equations. I have some code for Clairault, Lagrange and Riccati
equations which seems to work. I'd like some advice on integrating it
into maxima.
ODE2() is translated to lisp and compiled into maxima. I have done the
same for the new code, but put each new routine in its own file.
- ODE2() is modified to add a hook to the new code
- new routine ode1_nonlinear() calls
- ode1_clairault()
- ode1_lagrange()
- ode1_riccati()
So the questions are:
- is this structure reasonable
- should the compiled lisp code go into
- src, or
- share/diffequations
- if I polish the code and add the test cases will it be considered for maxima
--- ode2.mac.orig 2003-03-25 23:28:30.000000000 +1100
+++ ode2.mac 2003-04-07 23:42:27.000000000 +1000
@@ -44,7 +44,7 @@
ODE1A(EQ,Y,X):=BLOCK([DE,DES], /* f, g, %q% */
IF DERIVDEGREE(DE: EXPAND(LHS(EQ)-RHS(EQ)),Y,X) # 1
THEN RETURN(FAILURE(MSG1,EQ)),
- IF LINEAR2(DE,'DIFF(Y,X)) = FALSE THEN RETURN(FAILURE(MSG2,EQ)),
+ IF LINEAR2(DE,'DIFF(Y,X)) = FALSE THEN GO(NONLINEAR),
DES: DESIMP(DE),
DE: SOLVE1(DES,'DIFF(Y,X)),
IF FTEST(SOLVELNR(DE)) THEN RETURN(%Q%),
@@ -55,7 +55,10 @@
IF FTEST(INTEGFACTOR(%G%,%F%)) THEN RETURN(EXACT(%Q%*%G%,%Q%*%F%)),
IF FTEST(SOLVEHOM(DE)) THEN RETURN(%Q%),
IF FTEST(SOLVEBERNOULLI(DE)) THEN RETURN(%Q%),
- IF FTEST(GENHOM(DE)) THEN RETURN(%Q%) ELSE RETURN(FALSE))$
+ IF FTEST(GENHOM(DE)) THEN RETURN(%Q%),
+ NONLINEAR,
+ IF FTEST(ode1_nonlinear(DE,Y,X)) THEN RETURN(%Q%)
+ ELSE RETURN(FALSE))$
DESIMP(EQ):=BLOCK([INFLAG:TRUE],
EQ: FACTOR(EQ),
--- /dev/null 2003-04-08 12:00:28.000000000 +1000
+++ ode1_nonlinear.mac 2003-04-07 23:43:19.000000000 +1000
@@ -0,0 +1,19 @@
+/* -*- Mode: MACSYMA; Package: MAXIMA -*- */
+
+/* Nonlinear first order ordinary differential equations
+ This routine is called from ODE1A in ode2.mac, which
+ does the required checks*/
+
+declare(method,special);
+
+ode1_nonlinear(de,y,x):=block(
+ [q],
+ q:ode1_clairault(de,y,x),
+ if (is(q#false)) then return(q),
+ q:ode1_lagrange(de,y,x),
+ if (is(q#false)) then return(q),
+ q:ode1_riccati(de,y,x),
+ if (is(q#false)) then return(q),
+ return(false))$
+
+
--- /dev/null 2003-04-08 12:00:38.000000000 +1000
+++ ode1_clairault.mac 2003-04-08 00:19:12.000000000 +1000
@@ -0,0 +1,100 @@
+/* Solve Clairault's equation f(x*y'-y)=g(y')
+
+ References:
+
+ Daniel Zwillinger, Handbook of Differential Equations, 3rd ed
+ Academic Press, (1997), pp 216-218
+
+ Note: Singular solutions may also exist
+
+Test case
+ (x^2-1)*'diff(y,x)^2 - 2*x*y*'diff(y,x) + y^2 -1 = 0
+Solution
+ (%C*x-y)^2 = -%C^2-1
+
+Copyright (C) 2003 David Billinghurst
+
+This library is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+This library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with this library; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+*/
+
+declare(method,special);
+
+/* Attempt to seperate expression exp into f(x)+g(y)
+ Returns [f(x),g(y)] is this is possible
+ false otherwise
+ Adapted from SEPERABLE() in ode2.mac */
+
+plusseperable(eq,x,y) := block(
+ [u,xpart:[],ypart:[],flag:false,inflag:true],
+ eq:expand(eq),
+ if atom(eq) or not(inpart(eq,0)="+") then eq: [eq],
+ for u in eq do
+ if freeof(x,u) then
+ ypart: cons(u,ypart)
+ else if freeof(y,u) then
+ xpart: cons(u,xpart)
+ else
+ return(flag: true),
+ if flag = true then return(false),
+ if xpart = [] then xpart: 0 else xpart: apply("+",xpart),
+ if ypart = [] then ypart: 0 else ypart: apply("+",ypart),
+ return([xpart,ypart])
+)$
+
+/* If we can write eqn as f(x*y'-y)=g(y') then the solution
+ is given implicitly by f(x*%c-y)=g(%c)
+
+ A singular solution may exist
+
+ */
+ode1_clairault(eq,y,x) := block(
+ [ans,sing_ans,%C,de,%f,%g,%p,%r,%t,sep],
+ /* substitute %p=y', %r=x*y'-y */
+ de: expand(lhs(eq)-rhs(eq)),
+ de: subst(%p,'diff(y,x),de),
+ de: ratsimp(subst((%r+y)/%p,x,de)),
+ if not(freeof(x,y,de)) then return(false),
+ if (sep:plusseperable(de,%r,%p))=false then return(false),
+ %f: sep[1], %g:-sep[2],
+ method: 'CLAIRAULT,
+ ans:subst(x*%C-y,%r,%f)=subst(%C,%p,%g),
+ %t:solve(ans,y),
+ if length(%t)#0 then ans:first(%t),
+ sing_ans:ode1_clairault_singular(eq,y,x),
+ if length(sing_ans)#0 then ans:cons(ans,sing_ans),
+ return(ans)
+)$
+
+/* Is there a singular solution to Clairault equation eq:f(x*p-y)-g(p)=0 ?
+
+ Differentiate eq wrt x to get 'diff(y,x,2)*eq2=0
+ Eliminate 'diff(y,x) between eq and eq2.
+*/
+
+ode1_clairault_singular(eq,y,x) := block(
+ [expr,eq2,p,t],
+ eq:lhs(eq)-rhs(eq),
+ expr:subst(y(x),y,eq),
+ expr:diff(expr,x),
+ expr:subst(y,y(x),expr),
+ expr:expand(expr),
+ eq2:ratsimp(ratcoeff(expr,'diff(y,x,2))),
+ if not(freeof('diff(y,x,2),eq2)) then return([]),
+ if not(is(ratsimp(expr-eq2*'diff(y,x,2)=0))) then return([]),
+ /* eq2 is ok. return parametric solution as too hard to simplify */
+ eq:subst(t,'diff(y,x),eq=0),
+ eq2:subst(t,'diff(y,x),eq2=0),
+ return([[eq,eq2]])
+)$
--- /dev/null 2003-04-08 12:00:49.000000000 +1000
+++ ode1_lagrange.mac 2003-04-07 23:44:03.000000000 +1000
@@ -0,0 +1,40 @@
+
+/* Solution of Lagrange equation y = x f(y') + g(y')
+ Ref: Zwillinger (1st ed), pp265-268
+
+Test cases
+ (1) y=2*x*'diff(y,x)-a*'diff(y,x)^3
+ 2 4 3 2 4 2 2 3 2
+ - a (27ay +(-16 x -144%Cax)y +64%Cx +128%C ax +64%C a ) = 0
+
+*/
+
+declare(method,special);
+
+ode1_lagrange(eq,y,x):= block(
+ [eqn,t,f,g,t2,s,p,q],
+
+ /* Express equation as y = x*f(p)+g(p), where p=y' */
+ eqn:subst(p,'diff(y,x),eq),
+ if freeof(y,eqn) then return(false),
+ eqn:solve(eqn,y),
+ if eqn=[] then return(false),
+ eqn: eqn[1],
+ if lhs(eqn)#y then return(false),
+ if not(freeof(y,rhs(eqn))) then return(false),
+ t:rhs(eqn),
+ f:ratcoeff(t,x),
+ g:t-f*x,
+ if not(freeof(x,y,f)) then return(false),
+ if not(freeof(x,y,g)) then return(false),
+ if f=p then return(false),
+
+ t2:'diff(x,p)=x*diff(f,p)/(p-f)+diff(g,p)/(p-f),
+ s:ode2(t2,x,p),
+ if (s=false) then return(false),
+
+ q:first(eliminate([y=x*f+g,s],[p])),
+ if not(freeof(p,q)) then return(false),
+ method:'lagrange,
+ return(q=0)
+)$
--- /dev/null 2003-04-08 12:01:00.000000000 +1000
+++ ode1_riccati.mac 2003-04-08 00:00:20.000000000 +1000
@@ -0,0 +1,65 @@
+/* Attempt to solve Riccati ode y' = f2(x)*y^2+f1(x)*y+f0(x)
+
+ References:
+
+ D Zwillinger, Handbook of Differential Equations, 3rd ed
+ Academic Press, (1997), pp 354-355
+
+ G M Murphy, Ordinary Differential Equations and Their
+ Solutions, Van Nostrand, 1960, pp 15-23
+
+ Substitute y = -z'/(z*r)
+
+ => f2*z''-(f2'+f1*f2)z'+f2^2*f0*z=0
+
+ Solve this second order linear ode for z. The solution
+ has form z=%k1*f+%k2*g, with two constants, but a first
+ order ode only has one constant %c. Take %k1=1 and %k2=%c
+
+ y = -z'/(z*r) = -(f'+%c*g')/((f+%c*g)*r)
+
+Copyright (C) 2003 David Billinghurst
+
+This library is free software; you can redistribute it and/or
+modify it under the terms of the GNU Lesser General Public
+License as published by the Free Software Foundation; either
+version 2.1 of the License, or (at your option) any later version.
+
+This library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+Lesser General Public License for more details.
+
+You should have received a copy of the GNU Lesser General Public
+License along with this library; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+*/
+
+declare(method,special);
+
+ode1_riccati(eq,y,x) := block(
+ [de,%a,f0,f1,f2,ans,z,%C,%K1,%K2],
+ de:expand(lhs(eq)-rhs(eq)),
+ %a:coeff(de,'diff(y,x),1),
+ if %a=0 then return(false),
+ de:expand(de/%a),
+ f2:-ratsimp(coeff(de,y,2)),
+ if not(freeof(y,f2)) then return(false),
+ if f2=0 then return(false),
+ f1:-ratsimp(coeff(de,y,1)),
+ if not(freeof(y,f1)) then return(false),
+ f0:-ratsimp(de-'diff(y,x)+f2*y^2+f1*y),
+ if not(freeof(y,f0)) then return(false),
+ if not(is(de-'diff(y,x)+f2*y^2+f1*y+f0=0)) then return(false),
+
+ /* try to solve the new differential equation */
+ de: f2*'diff(z,x,2)-(diff(f2,x)+f1*f2)*'diff(z,x)+f2^2*f0*z=0,
+ if (ans:ode2(de,z,x))=false then return(false),
+
+ /* give up on parametric solutions for z */
+ if lhs(ans)#z then return(false),
+
+ method:'RICCATI,
+ ans:rhs(ans), ans:subst(1,%K1,ans), ans:subst(%C,%K2,ans),
+ return(y=ratsimp(-diff(ans,x)/(ans*f2)))
+)$