From: David B. <bil...@us...> - 2004-01-27 06:12:37
|
Update of /cvsroot/maxima/maxima/share/contrib/diffequations In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv11561 Added Files: ode1_riccati.mac Log Message: Initial version --- NEW FILE: ode1_riccati.mac --- /* ode1_riccati.mac 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 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_riccati,001,'version)$ ode1_riccati(eq,y,x) := block( [de,%a,f0,f1,f2,ans], 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:-expand(ratsimp(coeff(de,y,1))), if not(freeof(y,f1)) then return(false), f0:-expand(ratsimp(de-'diff(y,x)+f2*y^2+f1*y)), if not(freeof(y,f0)) then return(false), if not(is(ratsimp(expand(de-'diff(y,x)+f2*y^2+f1*y+f0))=0)) then return(false), ode_disp(" is Ricatti equation"), /* Following Murphy, (3-1, p15) see if the equation has the form of the original equation studied by Riccati y' + b*y^2 = c*x^m => f1 = 0 b:-f2(x) is constant f0 = c*x^m */ ans: block( [b,c,m], if ( f1#0 ) then return(false), if (not freeof(x,f2) ) then return(false), b:-f2, m:hipow(f0,x), c:coeff(f0,x,m), if ( ratsimp(f0-c*x^m)#0 ) then return(false), ode_disp(" equation is original Riccati equation"), ode1_riccati_original(b,c,m,y,x) ), if ( ans#false ) then ( method:'RICCATI, return(ans) ), /* Perhaps it has the special form Murphy (3-3, p21-22) x*y' -a*y + b*y^2 = c*x^n => a: f1(x)*x constant b:-f2(x)*x constant f0 = c*x^(n-1) */ ans: block( [a,b,c,m,n], if ( not freeof(x,a:ratsimp( f1*x)) ) then return(false), if ( not freeof(x,b:ratsimp(-f2*x)) ) then return(false), m:hipow(f0,x), c:coeff(f0,x,m), /* May want to check c and m */ if ( is(ratsimp(f0-c*x^m)#0) ) then return(false), n:m+1, ode_disp(" equation is special Riccati equation"), ode1_riccati_special(a,b,c,n,y,x) ), if ( ans#false ) then ( method:'RICCATI, return(ans) ), /* The equation doesn't have a special form, so it is a general Riccati equation. */ ans:ode1_riccati_general(f0,f1,f2,y,x), if ( ans#false ) then ( method:'RICCATI, return(ans) ), /* Default return value */ false )$ /* Solve the original Riccati equation y' + b*y^2 = c*x^m Murphy (3-2, p20-21) */ ode1_riccati_original(b,c,m,y,x) := block( [ans,k,w,s,i], ode_disp(" -> In ode1_riccati_original"), ode_disp2(" b: ",b), ode_disp2(" c: ",c), ode_disp2(" m: ",m), /* Solve m*(2*k+1)+4*k=0 => k= m/(2*m+4) If k is an integer then the equation is integrable in finite terms. The solution is then found using the transformation y = w/x, giving x*w'(x)-w+b*w^2=c*x^(m+2) which is the special Riccati equation with a=1 and n=m+2 */ if ( is(equal(m,-2)) ) then ( /* Case m=-2 is special. Do it first to avoid division by zero below */ ode1_riccati_original_2(b,c,m,y,x) ) else if ( integerp(k:m/(2*m+4)) ) then ( ode_disp(" Equation is integrable in finite terms"), ode_disp(" Transforming using y=w/x and calling ode1_riccati_special"), ans:ode1_riccati_special(1,b,c,m+2,w,x), if ( ans=false ) then error("Error in ode1_riccati_original"), ode_disp2(" Solution to transformed equation is ",ans), return(y=ratsimp(rhs(ans)/x)) ) else ( /* Transform to linear second order equation. Solution involves Bessel functions */ ode1_riccati_original_3(b,c,m,y,x) ) )$ /* Solve the original Riccati equation y' + b*y^2 = c*x^m for m=-2 - Transform to second order linear ode, Murphy (3-2c, p20-21) - Solve using Murphy A3-250 */ ode1_riccati_original_2(b,c,m,y,x) := block( ode_disp(" Special Riccati equation with m=-2"), ode_disp(" Equation is not integrable in finite terms"), ode_disp(" Solution not yet implemented"), false )$ /* Solve the original Riccati equation y' + b*y^2 = c*x^m for m#2 - Transform to second order linear ode, Murphy (3-2c, p20-21) - Solve using Murphy A3-41 */ ode1_riccati_original_3(b,c,m,y,x) := block( [p:m+2,n:1/(m+2),%c], ode_disp(" -> In ode1_riccati_original_3"), ode_disp(" Special Riccati equation with m#2"), ode_disp2(" b: ",b), ode_disp2(" c: ",c), ode_disp2(" m: ",m), ode_disp2(" p: ",p), ode_disp2(" n: ",n), ode_disp(" Not implemented"), return(false), /* Partial implementation: - answers not quite right - what if b<0 and/or c<0 */ /* Let b*y*u(x)=u'(x) so that ode becomes u''(x) + b*c*x^m*u(x) = 0 Solution expressed in terms of Bessel functions of order n */ if integerp(n) then ( u: sqrt(x)*(bessel_j(n,2*sqrt(b)*sqrt(c)*x^(p/2)/p) + %c*bessel_y(n,2*sqrt(b)*sqrt(c)*x^(p/2)/p) ) ) else ( u: sqrt(x)*(bessel_j(n,2*sqrt(b)*sqrt(c)*x^(p/2)/p) + %c*bessel_j(-n,2*sqrt(b)*sqrt(c)*x^(p/2)/p) ) ), return(y=ratsimp(diff(u,x)/(b*u))) )$ /* Solve the special Riccati equation x*y' -a*y + b*y^2 = c*x^n Murphy (3-3, p21-22) */ ode1_riccati_special(a,b,c,n,y,x) := block( [k,s], ode_disp(" -> In ode1_riccati_special"), ode_disp2(" a: ",a), ode_disp2(" b: ",b), ode_disp2(" c: ",c), ode_disp2(" n: ",n), /* Certain cases are integrable. */ /* Case (a.i). n=2*a Equation can be made exact using integrating factor x^(a-1) and integrated */ if ( is(equal(n,2*a)) ) then ( ode_disp(" Case (a.i)"), return(ode1_riccati_special_i(a,b,c,n,y,x)) ) /* Case (a.ii) (n-2*a)/(2*n) a positive integer */ else if ( integerp(k:(n-2*a)/(2*n)) and k>0 ) then ( ode_disp2(" Case (a.ii) with k = ",k), if oddp(k) then s:ode1_riccati_special_i(n/2,c,b,n,y,x) else s:ode1_riccati_special_i(n/2,b,c,n,y,x), if s#false then ( s:rhs(s), for i:(k-1) thru 1 step -1 do ( if oddp(i) then s:(a+i*n)/c+x^n/s else s:(a+i*n)/b+x^n/s ), return(y=a/b+x^n/s) ) ) /* Case (a.iii) (n+2*a)/(2*n) a positive integer */ else if ( integerp(k:(n+2*a)/(2*n)) and k>0 ) then ( ode_disp2(" Case (a.iii) with k = ",k), if oddp(k) then s:ode1_riccati_special_i(n/2,c,b,n,y,x) else s:ode1_riccati_special_i(n/2,b,c,n,y,x), if s#false then ( s:rhs(s), for i:(k-1) thru 1 step -1 do ( if oddp(i) then s:(i*n-a)/c+x^n/s else s:(i*n-a)/b+x^n/s ), return(y=x^n/s) ) ), /* Default return value */ ode_disp(" Equation is not integrable in finite terms"), ode_disp(" Solution can be expressed using Bessel functions"), ode_disp(" Not yet implemented"), false )$ /* Solve the special Riccati equation x*y' -a*y + b*y^2 = c*x^n for the case n=2*a. Murphy (3-3, p21-22). Note: Signs changed from Murphy in cases a.i.1 and a.i.3. */ ode1_riccati_special_i(a,b,c,n,y,x) := block( [%c,signb,signc], ode_disp(" -> In ode1_riccati_special_i"), ode_disp2(" a: ",a), ode_disp2(" b: ",b), ode_disp2(" c: ",c), ode_disp2(" n: ",n), if not(equal(n,2*a)) then error("ode1_riccati_special_i: n#2*a"), /* Case (a.i). n=2*a Equation can be made exact using integrating factor x^(a-1) and integrated to give solution(s) below. */ signb:asksign(b), signc:asksign(c), /* b*c > 0 */ if ( signb='pos and signc='pos ) then ( /* Murphy has the sign of the solution wrong */ ode_disp(" Case (a.i.1) b*c>0, b>0 and c>0"), return(y=sqrt(c/b)*x^a*tanh(sqrt(b*c)*x^a/a+%c)) ) else if ( signb='neg and signc='neg ) then ( ode_disp(" Case (a.i.2) b*c>0, b<0 and c<0"), return(y=sqrt(c/b)*x^a*tanh(%c-sqrt(b*c)*x^a/a)) ) /* b*c < 0 */ else if ( signb='pos and signc='neg ) then ( ode_disp(" Case (a.i.3) b*c<0, b>0 and c<0"), return(y=sqrt(-c/b)*x^a*tan(%c-sqrt(-b*c)*x^a/a)) ) else if ( signb='neg and signc='pos ) then ( ode_disp(" Case (a.i.4) b*c<0, b<0 and c>0"), /* Murphy has the sign of the solution wrong */ return(y=sqrt(-c/b)*x^a*tan(sqrt(-b*c)*x^a/a+%c)) ) /* b and c are non-zero constants, so this is an error */ else ( error("ode_riccati_special_i: Impossible case has just happened"), return(false) ) )$ /* The equation doesn't have a special form, so it is a generalized Riccati equation. Try transforming it to a linear second order ode. Substitute y = -z'/(z*f2) => 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. Without loss of generality take %k1=1 and %k2=%c y = -z'/(z*f2) = -(f'+%c*g')/((f+%c*g)*f2) */ ode1_riccati_general(f0,f1,f2,y,x) := block( [de,z,ans,%c,%k1,%k2], ode_disp(" Transforming to 2nd order ode"), de: f2*'diff(z,x,2)-(diff(f2,x)+f1*f2)*'diff(z,x)+f2^2*f0*z=0, if get('contrib_ode,'verbose) then disp(de), ans:contrib_ode(de,z,x), ode_disp(" with solution"), if get('contrib_ode,'verbose) then disp(ans), if is(ans=false) then ( ode_disp(" Cannot solve 2nd order ode"), false ) else if (lhs(ans[1])#z) then ( /* give up on parametric solutions for z */ ode_disp(" 2nd order ode has parametric solution"), ode_disp(" giving up"), false ) else ( ode_disp(" solution found"), method:'RICCATI, ans:rhs(ans[1]), ans:subst(1,%K1,ans), ans:subst(%C,%K2,ans), return(y=ratsimp(-diff(ans,x)/(ans*f2))) ) )$ |