Update of /cvsroot/maxima/maxima/share/contrib/diffequations In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv23133 Added Files: contrib_ode.mac contrib_ode.usg ode1_clairault.mac ode1_factor.mac ode1_lagrange.mac ode1_lie.mac ode1_nonlinear.mac Log Message: Initial versions --- NEW FILE: contrib_ode.mac --- /* contrib_ode.mac Driver for contributed ODE routines Copyright (C) 2004 David Billinghurst This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program 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 General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ if get('ode1_nonlinear,'version)=false then load('ode1_nonlinear)$ contrib_ode(eqn,y,x) := block( [soln,degree:derivdegree(eqn,y,x)], ode_disp("-> ode_contrib"), /* Try ode2 if first or second order*/ if (degree < 1) then ( ode_disp(" Degree of equation less than 1"), false ) else if (degree > 2) then ( ode_disp(" Degree of equation greater then 2"), false ) else if (degree=1) then ( ode_disp(" First order equation"), ode_disp(" -> ode2"), soln:ode2(eqn,y,x), if (soln#false) then ( ode_disp(" Successful"), return([soln]) ) else ( return(ode1_nonlinear(eqn,y,x)) ) ) else if (degree=2) then ( ode_disp(" Second order equation"), ode_disp(" -> ode2"), soln:ode2(eqn,y,x), if (soln#false) then ( ode_disp(" Successful"), [soln] ) else false ) else ( error("contrib_ode: This case cannot occur"), false ) )$ /* Check the solution of an ode. */ ode_check(e_,a_) := block( ratsimp(ev(lhs(e_)-rhs(e_),a_,diff)) )$ ode_disp(msg) := block( if get('contrib_ode,'verbose) then print(msg) )$ ode_disp2(msg,e) := block( if get('contrib_ode,'verbose) then print(msg,e) )$ --- NEW FILE: contrib_ode.usg --- -*- mode: Text; -*- MAXIMA's ordinary differential equation (ODE) solver ODE2 solves elementary linear ODEs of first and second order. The function contrib_ode extends ODE2 with additional methods. The code may be integrated into MAXIMA at some stage. The calling convention for contrib_ode is identical to ODE2. It takes three arguments: an ODE (only the left hand side need be given if the right hand side is 0), the dependent variable, and the independent variable. When successful, it returns a list of solutions. Each solution can be: an explicit solution for the dependent variable; an implicit solution for the dependent variable; or a parametric solution in terms of variable %t. %C is used to represent the constant in the case of first order equations, and %K1 and %K2 the constants for second order equations. If contrib_ode cannot obtain a solution for whatever reason, it returns FALSE, after perhaps printing out an error message. (C1) eqn:x*'diff(y,x)^2-(1+x*y)*'diff(y,x)+y=0; dy 2 dy (D1) x (--) - (x y + 1) -- + y = 0 dx dx (C2) contrib_ode(eqn,y,x); x (D2) [y = LOG(x) + %C, y = %C %E ] (C3) method; (D3) FACTOR Nonlinear odes can have singular solutions without constants of integration. (C4) eqn:'diff(y,x)^2+x*'diff(y,x)-y=0; dy 2 dy (D4) (--) + x -- - y = 0 dx dx (C5) contrib_ode(eqn,y,x); 2 2 x (D5) [y = %C x + %C , y = - --] 4 (C6) method; (D6) clairault The following ode has two parametric solutions in terms of the dummy variable %t. In this case the parametric solutions can be manipulated to give explicit solutions. (C7) eqn:'diff(y,x)=(x+y)^2; dy 2 (D7) -- = (y + x) dx (C8) contrib_ode(eqn,y,x); (D8) [[x = %C - ATAN(SQRT(%T)), y = - x - SQRT(%T)], [x = ATAN(SQRT(%T)) + %C, y = SQRT(%T) - x]] (C9) method; (D9) lagrange For first order ODEs contrib_ode calls ode2. It then tries the following methods: factorization, Clairault, Lagrange, Riccati and Lie symmetry methods. For second order ODEs contrib_ode calls ode2. No other methods are implemented yet. Extensive debugging traces and messages are displayed if the command put('contrib_ode,true,'verbose) is executed. TO DO ===== These routines are work in progress. I still need to: * Extend the FACTOR method ode1_factor to work for multiple roots. * Extend the FACTOR method ode1_factor to attempt to solve higher order factors. At present it only attemps to solve linear factors. * Fix the LAGRANGE routine ode1_lagrange to prefer real roots over complex roots. * Add additional methods for Riccati equations. * Work out how to return partial solutions, such as those obtained for Riccati equations. * Add routines for Abel equations. * Work on the Lie symmetry group routine ode1_lie. There are quite a few problems with it: some parts are unimplemented; some test cases seem to run forever; other test cases crash; yet others return very complex "solutions". I wonder if it really ready for release yet. * Add more test cases. TEST CASES ========== The routines have been tested on a few hundred test cases from Murphy, Kamke and Zwillinger. These are included in the tests subdirectory. * The Clairault routine ode1_clairault finds all known solutions, including singular soultions. (This statement is asking for trouble). * The other routines often return a single solution when multiple solutions exist. * Some of the "solutions" from ode1_lie are overly complex and impossible to check. * There are some crashes. References ========== [1] E Kamke, Differentialgleichungen Losungsmethoden und Losungen, Vol 1, Geest & Portig, Leipzig, 1961 [2] G M Murphy, Ordinary Differential Equations and Their Solutions, Van Nostrand, New York, 1960 [3] D Zwillinger, Handbook of Differential Equations, 3rd edition, Academic Press, 1998 --- NEW FILE: ode1_clairault.mac --- /* ode1_clairault.mac 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 Copyright (C) 2004 David Billinghurst This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program 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 General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ declare(method,special); put('ode1_clairault,001,'version)$ /* Attempt to separate expression exp into f(x)+g(y) Returns [f(x),g(y)] is this is possible false otherwise Adapted from SEPARABLE() in ode2.mac */ plusseparable(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,%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 ( ode_disp2(" Expression not free of x and y: ",de), return(false) ), sep:plusseparable(de,%r,%p), if is(sep=false) then /* Perhaps de is const+f(%p)*%r */ sep:plusseparable(ratsimp(de/%r),%r,%p), if is(sep#false) then ( %f: sep[1], %g:-sep[2] ) else ( ode_disp2(" Expression free of x and y but can't write as f(%r)=g(%p): ",de), /* Don't give up yet */ _t:solve(de,%r), if _t=[] then return(false), ode_disp2(" Solving for %r: ",_t), %f:%r, %g:rhs(_t[1]), if not(freeof(%r,%g)) then return(false) ), /* Next line added to avoid testsuite failures 2003-01-02 Guessed the fix and should analyse it further */ if ( is(%f=0) ) then return(false), ode_disp2(" %f: ",%f), ode_disp2(" %g: ",%g), method: 'clairault, ans:subst(x*%C-y,%r,%f)=subst(%C,%p,%g), /* Try and solve for y */ _t:solve(ans,y), if length(_t)=0 then ans:[ans] else ans:_t, ans:append(ans,ode1_clairault_singular(eq,y,x)), 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 If possible, eliminate 'diff(y,x) between eq and eq2. else return parametric solution in variable %t */ ode1_clairault_singular(eq,y,x) := block( [expr,eq2,%t,ans], 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([]), /* Try and eliminate t */ eq:subst(%t,'diff(y,x),eq), eq2:subst(%t,'diff(y,x),eq2), ans:eliminate([eq,eq2],[%t]), if ( not(freeof(%t,ans)) or (ans=1) ) then return([[eq=0,eq2=0]]), return(solve(ans,y)) )$ --- NEW FILE: ode1_factor.mac --- /* Solution of first order ODEs by factoring References: Daniel Zwillinger, Handbook of Differential Equations, 3rd ed Academic Press, (1997), pp 265-266 Copyright (C) 2004 David Billinghurst This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program 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 General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ declare(method,special); put('ode1_factor,001,'version)$ ode1_factor(eq,y,x):= block( [de,%p,ans:[],s,inflag:true,powflag:true,u,factors:[]], ode_disp(" in ode1_factor"), /* substitute %p=y' */ de: expand(lhs(eq)-rhs(eq)), de: subst(%p,'diff(y,x),de), /* At present restricted to simple cases such as p^2+5p+6=0 What about repeated factors? */ /* if ( not(freeof(x,y,de)) ) then return(false), */ de:factor(de), /* fail if the ode wasn't simplified by factoring */ if ( not(inpart(de,0)="*") or is(length(de)<2) ) then ( ode_disp2(" cannot factor",de), return(false) ), /* Count the factors that are DEs. There may be constant terms Some test cases failed when this was combined with the loop below */ for u in de do ( if not(freeof(%p,u)) then ( factors:endcons(u,factors) ) ), if length(factors)<2 then ( ode_disp2(" cannot factor",de), return(false) ), ode_disp2(" and after factoring is ",de), for u in de do ( /* Do not continue for higher order terms (including repeated factors) */ if hipow(expand(u),%p)>1 then ( ode_disp2(" unsuitable term ",u), powflag:false ) ), if powflag=false then return(false), /* For each factor, try and solve ODE */ for u in factors do ( ode_disp2(" Factor ",u), if not(freeof(%p,u)) then ( s:ode2(subst('diff(y,x),%p,u)=0,y,x), ode_disp2(" has solution",s), if s#false then ans:endcons(s,ans) ) ), if ans#[] then ( method:'factor, ans ) else false )$ --- NEW FILE: ode1_lagrange.mac --- /* ode1_lagrange.mac Solution of Lagrange equation y = x f(y') + g(y') References: Daniel Zwillinger, Handbook of Differential Equations, 3rd ed Academic Press, (1997), pp 328-331 Copyright (C) 2004 David Billinghurst This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program 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 General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ declare(method,special); put('ode1_lagrange,001,'version)$ ode1_lagrange(eq,y,x):= block( [eqn,ans1,ans:[],e,t,f,g,t2,s,p,q,%t], ode_disp("In ode1_lagrange"), /* 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), /* There may be multiple solutions */ for e in eqn do ( ans1:block( if lhs(e)#y then return(false), if not(freeof(y,rhs(e))) then return(false), t:rhs(e), 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), /* Transform ode */ t2:'diff(x,p)=x*diff(f,p)/(p-f)+diff(g,p)/(p-f), ode_disp2("ode1_lagrange: Transform successful and new equation is ",t2), s:ode2(t2,x,p), if (s=false) then return(false), ode_disp("Transformed equation solved"), method:'lagrange, t:y=x*f+g, /* Simple test to prevent failures */ ode_disp("Need to eliminate p from"), ode_disp2(" Equation t: ",t), ode_disp2(" Solution s: ",s), /* Eliminate is only for polynomials. Punt that it returns 1 when called inappropriately */ q:first(eliminate([t,s],[p])), ode_disp(" and after eliminating p the result is"), ode_disp2(" q: ",q), if ( freeof(p,q) and (q#1) ) then ( q=0 ) else ( /* Return the parametric equation in %t=p */ [subst(%t,p,s),subst(%t,p,t)] ) ), if ans1#false then ans:endcons(ans1,ans) ), if ans=[] then false else ans )$ --- NEW FILE: ode1_lie.mac --- /* ode1_lie.mac Solution of first order ODEs using Lie symmetry methods References: [1] E. S. Cheb-Terrab, A. D. Roche, Symmetries and First Order ODE Patterns [2] E. S. Cheb-Terrab, T. Koloknikov, First Order ODEs, Symmetries and Linear Transformations Copyright (C) 2004 David Billinghurst This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program 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 General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ put('ode1_lie,001,'version)$ /* Solve first order ODE y' = Phi(x,y) using Lie's method of symmetries */ ode1_lie(phi,y,x) := block( [sym,xi,eta,mu,_r], ode_disp("In ode1_lie"), /* Check that phi is free of dy/dx */ if not(freeof('diff(y,x),phi)) then ( ode_disp("In ode1_lie: phi not free of dy/dx"), return(false) ), /* Look for symmetry generators [xi,eta] */ sym:ode1_a(phi,y,x), if (sym=false) then ( ode_disp("In ode1_lie: ode1_a returned false"), return(false) ), xi:sym[1], eta:sym[2], ode_disp2("In ode1_lie: xi = ",xi), ode_disp2("In ode1_lie: eta = ",eta), /* Check symmetry generators */ r:symtest(phi,xi,eta,y,x), ode_disp2("In ode1_lie: symtest = ",r), /* Derive integrating factor from [xi,eta] */ mu:lie_integrating_factor(phi,xi,eta), /* Solve ODE using itegrating factor */ lie_exact(phi,mu,y,x) )$ /* Determine if ODE y' = Phi(x,y) has linear symmetry of form xi = F(x), eta = P(x)*y + Q(x) If so, we can transform into an ODE with symmetry xi = F(x), eta = H(x) for some F,H using lie_FxHx(). We then transform the symmetry back to the original coordinate system. If x=t and y=y(u,t) then F (d/dt) + H (d/du) = F (d/dx + (dy/dt) d/dy) + H (dy/du) d/dy = F d/dx + (F (dy/dt) + H (dy/du)) d/dy so that xi = F eta = F(dy/dt) + H(dy/du) Return [xi,eta] if successful false otherwise */ ode1_a(phi,y,x) := block( [ %A, %Ay, %Ayy, %Ayyy, %phi_yy, %phi_yyy ], ode_disp("In ode1_a ..."), ode_disp2(" phi = ",phi), /* Using [2], eq (13),(14) */ /* Confirm that diff(phi,y,3)#0 */ %phi_yy: ratsimp(diff(phi,y,2)), %phi_yyy: ratsimp(diff(%phi_yy,y)), if is(%phi_yyy=0) then ( /* %phi_yyy=0 - Not covered in [2] */ ode_disp(" %phi_yyy = 0"), return(lie_FxHx(phi,y,x)) ), %A: ratsimp(%phi_yy/%phi_yyy), %Ay: ratsimp(diff(%A,y)), %Ax: ratsimp(diff(%A,x)), if is(%Ay=0) then /* Case 1: diff(A,y)=0 Change variables using y=A*u (15) */ block( ode_disp(" ..... Case 1 true"), ode_disp(" Change variables using y=A.u"), ode_disp2(" where A is ",%A), /* Here neweq is RHS only */ neweq: ratsimp((subst(%A*u,y,phi)-%Ax*u)/%A), ode_disp2(" transformed equation is ",neweq), /* Symmetry generators of u equation */ ode_disp(" Finding symmetries of transformed equation"), gen:radcan(lie_FxHx(neweq,u,x)), ode_disp2(" [xi,eta] = ",gen), if (gen=false) then false else block( /* transform [F,H] to [xi,eta] */ [_F:gen[1],_H:gen[2],dydt:%Ax*y/%A,dudt:%A], [_F, ratsimp(_F*dydt+_H*dudt)] ) ) else if is((%Ayy:ratsimp(diff(%Ay,y)))=0) then /* Case 2: diff(A,y,2)=0 and diff(A,y)#0 Change variables using u=ln(A) (17) */ block( [sub], ode_disp(" ..... Case 2 true"), ode_disp2(" A is ",%A), ode_disp2(" ln(A) is ",log(%A)), sub: radcan(solve(u=log(%A),y)[1]), /* Robust enough? */ ode_disp2(" ",sub), neweq: ratsimp(ev(phi,sub)-diff(rhs(sub),x))/diff(rhs(sub),u), neweq:radcan(neweq), ode_disp2(" neweq is ",neweq), /* Symmetry generators of u equation */ ode_disp(" Finding symmetries of transformed equation"), gen:lie_FxHx(neweq,u,x), if (gen=false) then false else block( /* transform [F,H] to [xi,eta] */ [_F:gen[1],_H:gen[2],dydt,dudt,_xi,_eta], _xi:_F, dydt:diff(rhs(sub),x), dudt:diff(rhs(sub),u), _eta: _F*dydt+_H*dudt, /* Now eliminate u from eta */ _eta:subst(log(%A),u,_eta), _eta:radcan(_eta), [ _xi, _eta ] ) ) else block( [%Axy,%eye,%Iy,%p,%px], %Axy: ratsimp(diff(%Ay,x)), %eye: ratsimp(%Axy/%Ayy), if linear2(y,%eye) then /* Case 3: diff(A,y,2)#0 and I=Axx/Ayy linear in y p=exp(integrate(Iy,x)) Substitute y = u/p */ block( [ _F, _H, _xi, _eta ], ode_disp(" ..... Case 3 true"), if get('contrib_ode,'verbose) then print(" A is ",%A), %Iy: diff(%eye,y), %Iy:ratsimp(%Iy), %Iy:trigsimp(%Iy), if not(freeof(y,%Iy)) then ( print("I_y not free of y in ode1_a - impossible"), print("I_y = ",%Iy), return(false) ), %p: radcan(exp(integrate(%Iy,x))), ode_disp2(" p is ",%p), %px:diff(%p,x), neweq: ratsimp(%p*subst(u/%p,y,phi)+%px*u/%p), ode_disp2(" neweq is ",neweq), /* Symmetry generators of u equation */ ode_disp(" Finding symmetries of transformed equation"), gen:lie_FxHx(neweq,u,x), ode_disp2(" [xi,eta] = ",gen), if (gen=false) then false else ( /* Now we need to transform xi and eta back to original equation. Have F(x).(d/dx) + H(x).(d/du) y = u/%p */ _F:radcan(gen[1]), _H:radcan(gen[2]), _xi:_F, ode_disp2(" _xi = ",_xi), _eta: ratsimp(( _H - _F*%px*y )/%p), ode_disp2(" _eta = ",_eta), [ _xi, _eta ] ) ) else /* No such symmetry exists */ false ) )$ /* Determine if ODE y' = Phi(x,y) has linear symmetry of form xi = F(x), eta = H(x) This uses the notation from [1] E. S. Cheb-Terrab, A. D. Roche, Symmetries and First Order ODE Patterns */ lie_FxHx(Phi,y,x) := block( [Phi_y, Phi_yy, Q, Q_y], ode_disp2("In lie_FxHx with Phi:",Phi), Phi_y: ratsimp(diff(Phi,y)), Phi_yy: ratsimp(diff(Phi_y,y)), ode_disp2("Phi_yy = ",Phi_yy), if (Phi_yy=0) then ( ode_disp("Phi_yy=0 in lie_FxHx. This case not handled yet"), return(false) ), Q:ratsimp(Phi_y/Phi_yy), /* (39) */ Q_y:ratsimp(diff(Q,y)), if (Q_y#0) then /* Case 1: Q_y # 0 */ block( [Phi_x,Upsilon,Upsilon_x,Upsilon_y,_t,_f,_h], Upsilon: ratsimp(diff(Q,x)/Q_y), /* (40) */ Upsilon_y: ratsimp(diff(Upsilon,y)), if is(Upsilon_y#0) then return(false), Phi_x: diff(Phi,x), Upsilon_x: diff(Upsilon,x), _t: ratsimp((Upsilon*Phi_y-Upsilon_x-Phi_x)/(Phi+Upsilon)), if is(ratsimp(diff(_t,y))#0) then return(false), _f: ratsimp(exp(integrate(_t,x))), /* (43) */ _h: ratsimp(-Upsilon*_f), [_f,_h] ) else /* Case 2: Q_y = 0 */ block( [_a,_b,_c,_f,_h], /* In this case Q is constant and Phi must have the form y' = Phi(x,y) = A(x) + B(x)*exp(y/Q) and F(x), H(x) are simple functions of A and B */ _c: bothcoef(expand(Phi),exp(y/Q)), _a: second(_c), _b: first(_c), if (_b=0) then ( print("Impossible case in lie_FxHx. This cannot happen."), error("Unknown error in lie_FxHx") ), _f: exp(-integrate(_a/Q,x))/_b, _h: _a*_f, [_f,_h] ) )$ /* Check if [ xi(x,y), eta(x,y) ] are symmetry generators of first order ode y'=phi(x,y) Returns 0 when the symmetry generators are OK. If the result is non-zero, then further simplification may reduce it to 0. */ symtest(phi,xi,eta,y,x) := ratsimp( diff(eta,x) + (diff(eta,y)-diff(xi,x))*phi - diff(xi,y)*phi^2 - xi*diff(phi,x) - eta*diff(phi,y) )$ /* Determine integrating factor of first order ode y'=phi(x,y) given symmetry generators [ xi(x,y), eta(x,y) ] */ lie_integrating_factor(phi,xi,eta) := 1/ratsimp(eta-xi*phi)$ /* Solve DE given integrating factor */ lie_exact(phi,mu,y,x) := block( [a,%c], a: integrate(mu*phi,x), method:'lie, intfactor: mu, ratsimp(-a+integrate(ratsimp(mu+diff(a,y)),y))=%c )$ /* from ode2.mac EXACT(M,N):=BLOCK([A,ynew,%c], INTFACTOR: SUBST(YOLD,YNEW,%q%), A: INTEGRATE(RATSIMP(M),X), METHOD: 'EXACT, RETURN(RATSIMP(A + INTEGRATE(RATSIMP(N-DIFF(A,Y)),Y)) = %C))$ */ --- NEW FILE: ode1_nonlinear.mac --- /* ode1_nonlinear.mac Driver routine for first order nonlinear ODE routines Copyright (C) 2004 David Billinghurst This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program 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 General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ declare(method,special); put('ode1_nonlinear,001,'version)$ if get('ode1_clairault,'version)=false then load("ode1_clairault.mac")$ if get('ode1_lagrange,'version)=false then load("ode1_lagrange.mac")$ if get('ode1_riccati,'version)=false then load("ode1_riccati.mac")$ if get('ode1_lie,'version)=false then load("ode1_lie.mac")$ if get('ode1_factor,'version)=false then load("ode1_factor.mac")$ ode1_nonlinear(de,y,x):=block( [q], ode_disp(" in ode1_nonlinear"), ode_disp(" -> ode1_factor"), q:ode1_factor(de,y,x), if (is(q#false)) then return(q), ode_disp(" -> ode1_clairault"), q:ode1_clairault(de,y,x), if (is(q#false)) then return(q), ode_disp(" -> ode1_lagrange"), q:ode1_lagrange(de,y,x), if (is(q#false)) then return(q), ode_disp(" -> ode1_riccati"), q:ode1_riccati(de,y,x), if (is(q#false)) then return([q]), ode_disp(" -> ode1_lie"), q:block( [phi], phi:solve(lhs(de)-rhs(de),'diff(y,x)), if (is(phi=[])) then false else ( phi:rhs(first(phi)), q:ode1_lie(phi,y,x) ) ), if (is(q#false)) then return([q]), return(false) )$ |