The commnand DEtools:-odeadvisor(ode,y(x)); gives the ode type names.
See https://maplesoft.com/support/help/maple/view.aspx?path=DEtools/odeadvisor for known names of ode type in Maple.
But we can also tell it to check if the ode is of specific type using
DEtools:-odeadvisor(ode,y(x),[name]);. For an example
ode:=x^2+3*x*diff(y(x),x)=y(x)^4+2*y(x); DEtools:-odeadvisor(ode,y(x),[Chini]) [_Chini]
If the name we have given is not the type of the ode, Maple will return [NONE].
For an example
ode:=x^2+3*x*diff(y(x),x)=y(x)^4+2*y(x); DEtools:-odeadvisor(ode,y(x),[separable]) [NONE]
To find what methods dsolve uses do
`dsolve/methods`[1] [quadrature, linear, Bernoulli, separable, inverse_linear, homogeneous, Chini, lin_sym, exact, Abel, pot_sym]
The above gives methods for first order ODEs. More are given by
`dsolve/methods`[1,'semiclass'] [Riccati, inverse_Riccati, equivalent_to_Abel, linearizable, linearizable_by_differentiation] `dsolve/methods`[1,'high_degree'] [WeierstrassP, WeierstrassPPrime, JacobiSN, linearizable_by_differentiation, missing, dAlembert, homogeneous_B, sym_implicit] `dsolve/methods`[1,"development"] [linearizable_by_differentiation, linearizable, con_sym, WeierstrassP, WeierstrassPPrime, equivalent_to_Abel, Abel_AIR, special, Riccati_symmetries] `dsolve/methods`[1,extra] [inverse_Riccati, Abel_AIL, `sym_pat/[F(x)*G(y),0]`, `sym_pat/[F(x),G(x)]`, `sym_pat/[F(x),G(y)]`, `sym_pat/[F(x)+G(y),0]`, `sym_pat/[F(x),G(x)*y+H(x)]`, sym_pat, exp_sym]
For example given a first order ode, we can ask dsolve to solve it using one of these methods as follows
restart; ode:=diff(y(x),x)=sqrt( sqrt( y(x) )* sin(x))/sqrt(y(x)) dsolve(ode,y(x),[`exact`]) 4/5*y(x)^(5/4) + Intat(-sqrt(sqrt(y(x))*sin(_a))/y(x)^(1/4), _a = x) + c__1 = 0
We can ask it to use symmetry with specific pattern as follows
restart; ode:=diff(y(x),x)=sqrt( sqrt( y(x) )* sin(x))/sqrt(y(x)) dsolve(ode,y(x),[`sym_pat/[F(x),G(x)*y+H(x)]`]) c__1 + 4/5*y(x)^(3/2)*sqrt(sin(x))/sqrt(sqrt(y(x))*sin(x)) - Int(sqrt(sin(x)), x) = 0 #or using short cut, same thing can be done as follows dsolve(ode,y(x),Lie) c__1 + 4/5*y(x)^(3/2)*sqrt(sin(x))/sqrt(sqrt(y(x))*sin(x)) - Int(sqrt(sin(x)), x) = 0 #another short cut for Lie is dsolve(ode,[`sym_pat`]) c__1 + 4/5*y(x)^(3/2)*sqrt(sin(x))/sqrt(sqrt(y(x))*sin(x)) - Int(sqrt(sin(x)), x) = 0
To find all methods, and for higher order ode, we first use indices as follows
indices(`dsolve/methods`) [high, linear_nonhomogeneous], [1], [1, high_degree], [3, linear_homogeneous], [2, "linear_homogeneous all methods"], [2, "linear_homogeneous other"], [2, linear_homogeneous], [2, "special_functions"], [3, linear_nonhomogeneous], [2, "linear_homogeneous in Normal Form"], [2, "development"], [2, "hypergeometric"], [1, extra], [high, linear_homogeneous], [high, nonlinear], [1, "special"], [1, "development"], [2, nonlinear], [high, "development"], [1, semiclass], [3, nonlinear], [2, linear_nonhomogeneous], [2, "linear_homogeneous as given"], [3, "development"]
Using the above, to find all methods for second order linear_homogeneous ode, the command is
`dsolve/methods`[2, "linear_homogeneous all methods"] [quadrature, const_coeffs, Euler, linear_1, `linear/missing_y`, Kovacic, special_functions, to_const_coeffs, exact_linear, sym_1, Mathieu, MeijerG, Heun, HeunG, HeunC, HeunB, HeunD, HeunT, mu_xy, equivalent_to_Bessel, to_Riccati, Bessel, elliptic, Legendre, Whittaker, Kummer, cylindrical, hypergeometric, hypergeom1, hypergeom2, Riemann, RNF, hypergeometricsols, rationalize_lode, with_periodic_functions] `dsolve/methods`[2, "linear_homogeneous other"] [exact_linear, sym_1, to_const_coeffs, mu_xy, equivalent_to_Bessel, to_Riccati, with_periodic_functions]
For example, given the ode \(y''+3 y'+ y=0\), we can now do
dsolve(ode,y(x)) y(x) = c__1*exp(1/2*(sqrt(5) - 3)*x) + c__2*exp(-1/2*(3 + sqrt(5))*x) dsolve(ode,y(x),[`const_coeffs`]) y(x) = c__1*exp(1/2*(sqrt(5) - 3)*x) + c__2*exp(-1/2*(3 + sqrt(5))*x) dsolve(ode,y(x),[`Kovacic`]) y(x) = c__1*exp(1/2*(sqrt(5) - 3)*x) + c__2*exp(-1/2*(3 + sqrt(5))*x)
Not all methods ofcourse will work, as it depends on the ode type.
This function below lists all methods
ind:=indices(`dsolve/methods`); for item in ind do cat("`dsolve/methods`",String(item)); eval(parse(%)) od;
Which gives
"`dsolve/methods`[high, linear_nonhomogeneous]" [quadrature, fully_exact_linear, linear_nonhomogeneous_[0,1], exact_linear_nonhomogeneous, linear, exp_reduce] "`dsolve/methods`[1]" [quadrature, linear, Bernoulli, separable, inverse_linear, homogeneous, Chini, lin_sym, exact, Abel, pot_sym] "`dsolve/methods`[1, high_degree]" [WeierstrassP, WeierstrassPPrime, JacobiSN, linearizable_by_differentiation, missing, dAlembert, homogeneous_B, sym_implicit] "`dsolve/methods`[3, linear_homogeneous]" [quadrature, const_coeffs, Euler, fully_exact_linear, to_const_coeffs, linear, exp_reduce, exact_linear, with_periodic_functions] "`dsolve/methods`[2, "linear_homogeneous all methods"]" [quadrature, const_coeffs, Euler, linear_1, linear/missing_y, Kovacic, special_functions, to_const_coeffs, exact_linear, sym_1, Mathieu, MeijerG, Heun, HeunG, HeunC, HeunB, HeunD, HeunT, mu_xy, equivalent_to_Bessel, to_Riccati, Bessel, elliptic, Legendre, Whittaker, Kummer, cylindrical, hypergeometric, hypergeom1, hypergeom2, Riemann, RNF, hypergeometricsols, rationalize_lode, with_periodic_functions] "`dsolve/methods`[2, "linear_homogeneous other"]" [exact_linear, sym_1, to_const_coeffs, mu_xy, equivalent_to_Bessel, to_Riccati, with_periodic_functions] "`dsolve/methods`[2, linear_homogeneous]" [linear_homogeneous] "`dsolve/methods`[2, "special_functions"]" [Bessel, elliptic, Legendre, Kummer, Whittaker, hypergeometric, Mathieu] "`dsolve/methods`[3, linear_nonhomogeneous]" [quadrature, fully_exact_linear, linear_nonhomogeneous_[0,1], exact_linear_nonhomogeneous, linear, exp_reduce] "`dsolve/methods`[2, "linear_homogeneous in Normal Form"]" [linear_1] "`dsolve/methods`[2, "development"]" [mu_xyp, mu_xyp2, mu_formal, mu_heuristic, exp_reduce, linear, Bessel2, Whittaker_old, nonlinear_homogeneous, exact_linear_nonhomogeneous, mu_y1, mu_x_y1, mu_y_y1, mu_poly_yn, exp_sym, sym_pat, sym_8] "`dsolve/methods`[2, "hypergeometric"]" [hypergeom1, hyper3] "`dsolve/methods`[1, extra]" [inverse_Riccati, Abel_AIL, sym_pat/[F(x)*G(y),0], sym_pat/[F(x),G(x)], sym_pat/[F(x),G(y)], sym_pat/[F(x)+G(y),0], sym_pat/[F(x),G(x)*y+H(x)], sym_pat, exp_sym] "`dsolve/methods`[high, linear_homogeneous]" [quadrature, const_coeffs, Euler, fully_exact_linear, to_const_coeffs, linear, exp_reduce, exact_linear, with_periodic_functions] "`dsolve/methods`[high, nonlinear]" [linearizable_by_differentiation, linearizable, reducible, exact_nonlinear, missing, mu_formal, lin_sym] "`dsolve/methods`[1, "special"]" [80, 81] "`dsolve/methods`[1, "development"]" [linearizable_by_differentiation, linearizable, con_sym, WeierstrassP, WeierstrassPPrime, equivalent_to_Abel, Abel_AIR, special, Riccati_symmetries] "`dsolve/methods`[2, nonlinear]" [Liouville, WeierstrassP, JacobiSN, linearizable, linearizable_by_differentiation, mu_xy_2, missing, mu_xyp2_dynamical_symmetries_fully_reducible, mu_xyp_singularcases, sym_1, exact_nonlinear, reducible, lin_sym, S-function, mu_xyp_generalcase, mu_xyp2_dynamical_symmetries_not_fully_reducible] "`dsolve/methods`[high, "development"]" [k25, RNF, mu_heuristic, MeijerG, nonlinear_homogeneous, mu_poly_yn, exp_sym] "`dsolve/methods`[1, semiclass]" [Riccati, inverse_Riccati, equivalent_to_Abel, linearizable, linearizable_by_differentiation] "`dsolve/methods`[3, nonlinear]" [linearizable_by_differentiation, linearizable, missing, exact_nonlinear, reducible, mu_formal, lin_sym] "`dsolve/methods`[2, linear_nonhomogeneous]" [quadrature, fully_exact_linear, linear_nonhomogeneous_[0,1], linear_nonhomogeneous_[0,F(x)], linear_nonhomogeneous] "`dsolve/methods`[2, "linear_homogeneous as given"]" [quadrature, const_coeffs, Euler, linear_1, linear/missing_y, Kovacic, RNF, special_functions, MeijerG, Heun, hypergeometricsols, rationalize_lode] "`dsolve/methods`[3, "development"]" [k25, RNF, mu_heuristic, linear_patterns, MeijerG, nonlinear_homogeneous, mu_y2, mu_poly_yn, exp_sym, pFq, 3F2, 2F2, 1F2, 0F2]
restart; ode:=diff(y(x),x)+y(x)^2*sin(x)-2*sin(x)/cos(x)^2 = 0; yp:=DETools:-particularsol(ode);
To step into the code, do
restart; ode:=diff(y(x),x)+y(x)^2*sin(x)-2*sin(x)/cos(x)^2 = 0; stopat(`DEtools/particularsol`); DETools:-particularsol(ode);
To print it do
print(`DEtools/particularsol`);
Use the output=basis option
ode:=diff(y(x),x$2)-x*diff(y(x),x)-x*y(x)=0; dsolve(ode,output=basis);
To solve
with \(y(0)=1,y'(0)=5\) do
eq1:= diff(y(x),x$2)-3*diff(y(x),x)+2*y(x)=10*exp(5*x); dsolve({eq1,y(0)=1,D(y)(0)=5},y(x)); Methods for second order ODEs: Trying to isolate the derivative d^2y/dx^2... Successful isolation of d^2y/dx^2 --- Trying classification methods --- trying a quadrature trying high order exact linear fully integrable trying differential order: 2; linear nonhomogeneous with symmetry [0,1] trying a double symmetry of the form [xi=0, eta=F(x)] <- double symmetry of the form [xi=0, eta=F(x)] successful ....
The above can also be written using D@@ notation, like this
eq:= (D@@2)(y)(x) - 3*D(y)(x) +2*y(x) = 10*exp(5*x); IC := y(0)=1,D(y)(0)=5; dsolve({eq,IC},y(x));
use odetest and check if it gives zero.
eq1:= diff(diff(y(x),x),x)-3*diff(y(x),x)+2*y(x)=10*exp(5*x); ans:=dsolve({eq1,IC},y(x)); odetest(ans,eq1); 0
Maple can classify the ODE.
eq1:= diff(y(x),x$2)-3*diff(y(x),x)+2*y(x)=10*exp(5*x); R0 := DEtools['odeadvisor'](eq1,y(x)); R0 := [[_2nd_order, _with_linear_symmetries]]
To get help on this type of ODE, do
DEtools['odeadvisor'](eq1,'help');
Use with(DEtools);
restart; eq1:= diff(y(x),x$2)-3*diff(y(x),x)+2*y(x)=10*exp(5*x); DEtools[DEplot](eq1,y(x),x=-2..5, [ [y(0)=0, D(y)(0)=0]], y=-3..3,linecolor=red);
To get a better plot, change the stepsize and independent variable range
restart; eq1:= diff(y(x),x$2)-3*diff(y(x),x)+2*y(x)=10*exp(5*x); DEtools[DEplot](eq1,y(x),x=-1..1,[[y(0)=0,D(y)(0)=0]],y=-3..3,stepsize=0.001,linecolor=red);
From: Robert Israel (israel@math.ubc.ca) Subject: Re: given precision in Maple Newsgroups: comp.soft-sys.math.maple Date: 2003-07-16 20:19:06 PST Set Digits:= n and all calculations from this point will be done with n digits. Mathematical functions will be correct to n digits as well (to the extent this is practical). If you want high-accuracy numerical ODE solutions, on the other hand, it's not so simple. I think the best way is using the taylorseries method. For example, consider the problem y' = y^2, y(1) = 1, where the exact solution y = 1/(2-x) has y(1.9) = 10. > Digits:= 30: sol:= dsolve({D(y)(x)=y(x)^2, y(1) = 1}, y(x), numeric, method=taylorseries, abserr=1e-25): sol(1.9); [x = 1.9, y(x) = 9.99999999999999999999999797691] > 10 - eval(y(x),%); -23 0.202309 10 The other methods (in particular the default rkf45) do not give results anywhere near this good.
How to obtain the ODE in canonical coordinates \(S(R)\) for an ode using Lie symmetry?
Here is a function I wrote with the help of Maple docs, that does that. For now, I am using first order ode. It takes the ode and the variables (coordinates) to use (\(R,S\)) is what normally used, and returns back the ODE in canonical coordinates. The ode should always be a quadrature for first order.
get_ODE_in_canonical:=proc(ode::`=`,y::symbol,x::symbol,S::symbol,R::symbol) local infinitesimals,tr,itr,ODE; infinitesimals := [DEtools:-symgen(ode)]; if nops(infinitesimals)=0 then RETURN(FAIL); fi; infinitesimals:=infinitesimals[1]; tr := DEtools:-canoni(infinitesimals,y(x),S(R)); itr := op(1,[solve(tr,{x,y(x)})]); ODE:= PDEtools:-dchange(itr,ode,[R,S(R)],simplify); ODE:= op(solve(ODE,{diff(S(R),R)})); RETURN(ODE); end proc:
Example
ODE := diff(y(x),x) = 1/x^2+y(x); get_ODE_in_canonical(ODE,y,x,S,R) #diff(S(R), R) = 1/(exp(R)*R^2) ODE := diff(y(x),x) = x^2-2*x*y(x)+y(x)^2; get_ODE_in_canonical(ODE,y,x,S,R) #diff(S(R), R) = 1
The following function is small change of the above. It accepts your choice of \(\xi ,\eta \) to use instead of the first one obtained by Maple. This makes it possible to experiment with different infinitesimals. Note that different infinitesimals will/can result in different canonical ODE, but the final solution for \(y(x)\) of course will always be the same after applying the transformation back to \(x,y\) coordinates.
restart; get_ODE_in_canonical_v2:=proc(ode::`=`,y::symbol,x::symbol,S::symbol,R::symbol,xi,eta) local infinitesimals:=[xi,eta],tr,itr,ODE; tr := DEtools:-canoni(infinitesimals,y(x),S(R)); itr := op(1,[solve(tr,{x,y(x)})]); ODE:= PDEtools:-dchange(itr,ode,[R,S(R)],simplify); ODE:= op(solve(ODE,{diff(S(R),R)})); RETURN(simplify(ODE,symbolic)); end proc:
And now
ODE:=diff(y(x),x) = x^2-2*x*y(x)+y(x)^2; get_ODE_in_canonical_v2(ODE,y,x,S,R,1,1) #diff(S(R), R) = 1/(R^2 - 1)
Here is an example.
ode:=diff(y(x),x)=-y(x)^2/(exp(x)-y(x)) DEtools:-symgen(ode,y(x),HINT=[c1*x+c2*y+c3,c4*x+c5*y+c6]); [_xi = 1, _eta = y]
Or
ode:=diff(y(x),x)=-y(x)^2/(exp(x)-y(x)) DEtools:-symgen(ode,y(x),HINT=[g(x),f(x)*y]); [_xi = 1, _eta = y]
So not use y(x) but only y in the hint.
The following function takes in a single ode and parses it to verify it is valid. It returns back 3 values. The dependent variable, the indepenent variable and the ode order.
#added 12/8/2022. interface(warnlevel=4); kernelopts('assertlevel'=2): parse_single_ode:=proc(ode::`=`)::symbol,symbol,integer; #parses single ode. returns back dep_var,indep_var,ode_order #throws exception when it detects parsing errors local func; local y,x; local dep_variables_found::list,item; local the_order; func:=PDEtools:-Library:-GetDepVars("not given",ode,onlydifferentiatedfunctions=true); if nops(func)=0 then error ("not differential equation ",ode); fi; func := func[1]; if nops(func)<>1 then error("Parsing error, dependent variable must contain one argument, found ", func); fi; y:=op(0,func); x:=op(1,func); #basic verification if not has(ode,y) then error ("Supplied ode ",ode," has no ",y); fi; if not has(ode,x) then error ("Supplied ode ",ode," has no ",x); fi; if not has(ode,func) then error ("Supplied ode ",ode," has no ",func); fi; the_order := PDEtools:-difforder(ode,x); #if the_order=0 then # error ("No derivative found in ",ode,". Input is not differential equation"); #fi; #note that the following call will also return y(x) if the input is not an ode #this will check that the dependent variable will show with SAME argument in the ode #i.e. if y(x) and y(t) show up in same ode, it will throw exception, which is what #we want. try dep_variables_found := PDEtools:-Library:-GetDepVars([y],ode); catch: error lastexception; end try; #now go over dep_variables_found and check the independent variable is same as x #i.e. ode can be y'(z)+y(z)=0 but function is y(x). for item in dep_variables_found do if not type(item,function) then error("Parsing error. Expected ",func," found ",item," in ode"); else if op(1,item) <> x then error("Parsing error. Expected ",func," found ",item," in ode"); fi; fi; od; #now go over all indents in ode and check that y only shows as y(x) and not as just y #as the PDEtools:-Library:-GetDepVars([_self:-y],ode) code above does not detect this. #i.e. it does not check y'(x)+y=0 if numelems(indets(ode,identical(y))) > 0 then error("Parsing error, Can not have ",y," with no argument inside ",ode); fi; return y,x,the_order; end proc:
To use do
y,x,the_order := parse_single_ode(diff(y(x),x)+y(x)=0);
An alternative to the above is to pass the dependent function itself as well as the ode. This is what I do myself in my ode solver. Like this
parse_single_ode:=proc(ode::`=`,func::function(name))::symbol,symbol,integer; #parses single ode. returns back dep_var,indep_var,ode_order #throws exception when it detects parsing errors local y,x; local dep_variables_found::list,item; local the_order; if nops(func)<>1 then error("Parsing error, dependent variable must contain one argument, found ", func); fi; y:=op(0,func); x:=op(1,func); #basic verification if not has(ode,y) then error ("Supplied ode ",ode," has no ",y); fi; if not has(ode,x) then error ("Supplied ode ",ode," has no ",x); fi; if not has(ode,func) then error ("Supplied ode ",ode," has no ",func); fi; the_order := PDEtools:-difforder(ode,x); #note that the following call will also return y(x) if the input is not an ode #this will check that the dependent variable will show with SAME argument in the ode #i.e. if y(x) and y(t) show up in same ode, it will throw exception, which is what #we want. try dep_variables_found := PDEtools:-Library:-GetDepVars([y],ode); #print("dep_variables_found=",dep_variables_found); catch: error lastexception; end try; #now go over dep_variables_found and check the independent variable is same as x #i.e. ode can be y'(z)+y(z)=0 but function is y(x). for item in dep_variables_found do if not type(item,function) then error("Parsing error. Expected ",func," found ",item," in ode"); else if op(1,item) <> x then error("Parsing error. Expected ",func," found ",item," in ode"); fi; fi; od; #now go over all indents in ode and check that y only shows as y(x) and not as just y #as the PDEtools:-Library:-GetDepVars([_self:-y],ode) code above does not detect this. #i.e. it does not check y'(x)+y=0 if numelems(indets(ode,identical(y))) > 0 then error("Parsing error, Can not have ",y," with no argument inside ",ode); fi; return y,x,the_order; end proc:
To use do the same as before, but need to add \(y(x)\) as second argument. Like this
y,x,the_order := parse_single_ode(diff(y(x),x)+y(x)=0, y(x));
This is very much the same as above, but this function returns True if the ode is syntactly valid, else false. Can be used just to check if the ode is valid before using it.
It takes in the ode and the function \(y(x)\) and returns either true or false.
#March 15, 2024 export is_valid_single_ode :=proc(ode_in::`=`,func::function(name))::truefalse; local x:=op(1,func),y:=op(0,func); local ode:=ode_in; local dep_variables_found::list; local item; if nops(func)<>1 then RETURN(false); fi; if not has(ode,diff) then ode:=convert(ode,diff); fi; if `or`(not has(ode,diff), not has(ode,x), not has(ode,func)) then RETURN(false); fi; try dep_variables_found := PDEtools:-Library:-GetDepVars([y],ode); catch: RETURN(false); end try; map( X-> `if`( `or`(not type(X,function), op(1,X) <> x) ,RETURN(false),NULL), dep_variables_found ); #check there is no y on its own. Should always be y(x) if nops(indets(ode,identical(y))) <> 0 then RETURN(false); fi; RETURN(true); end proc:
Example usages
is_valid_single_ode(diff(y(x),x)+sin(x)=0, y(x) ) true
To test, do
L:=[[3*(D@@2)(y)(x)+diff(y(x),x)+1=sin(x),y(x),true], [diff(y(x),x)^2=z,y(x),true], [diff(y(x),x)^2=z,y(x),true], [diff(y(x),x)^2=y(x)^2,y(x),true], [diff(y(x),x)=y(z)^2,y(x),false], [1/diff(y(x),x)=y(x)^2,y(x),true], [y(x)=1,y(x),false], [diff(y(x),z)=1,y(x),false], [D(y)(x)=1,y(x),true], [(D@@2)(y)(x)+diff(y(x),x)=sin(x),y(x),true], [diff(y(x),x)+diff(y(z),z)=1,y(x),false], [diff(y(x),x)+diff(g(z),z)=1,y(x),true] ]: map(X->evalb(is_valid_single_ode(X[1],X[2])=X[3]),L); if andmap(x->evalb(x=true),%) then print("all tests passed"); else print("WARNING, not all tests passed"); fi; #gives [true, true, true, true, true, true, true, true, true, true, true, true] "all tests passed"
Suppose to want to check if an ode has \(y'(x)\) in it. We can not write has(ode, diff(y(x),x))
and see if this gives true or not, because this will also match \(y''\) and \(y'''\) and any higher
order.
One way is to first convert the ode to D form and then use has on D(y)(x). This will not
match \(y''\) anymore which is what we wanted, because \(y''\) becomes (D@@2)(y)(x) and so the check
for first order diff will work as expected.
ode:=diff(y(x),x$2)+3*y(x)=0; ode:=convert(ode,D); has(ode,(D)(y)(x)) #false
Another example
ode:=diff(y(x),x$2)+1/diff(y(x),x)=0; ode:=convert(ode,D); has(ode,(D)(y)(x)) #true
The same thing if we wanted to check if the ode has \(y''\) or not.
ode:=diff(y(x),x$3)+diff(y(x),x)=0; ode:=convert(ode,D); has(ode,(D@@2)(y)(x)) #false
ode:=diff(y(x),x$3)+ x*diff(y(x),x)+sin(x)=0; if has(DEtools:-odeadvisor(ode,y(x),['linear']),_linear) then print("linear"); else print("not linear ode"); fi
The above prints "linear"
ode:=y(x)*diff(y(x),x$3)+ x*diff(y(x),x)+sin(x)=0; if has(DEtools:-odeadvisor(ode,y(x),['linear']),_linear) then print("linear"); else print("not linear ode"); fi
The above prints "not linear ode". This works for any ode order.
ode:=diff(y(x),x$3)+ x*diff(y(x),x)+sin(x)=0; PDEtools:-difforder(ode)
Gives \(3\).
Given a linear ode of the form \(A y''+B y' + C y= f(x)\) how to find \(A,B,C,f(x)\) ?
ode:=diff(y(x),x$3)+ x*diff(y(x),x)+99*y(x)=sin(x); L:=DEtools:-convertAlg(ode,y(x)); #this only works on linear ode's
Gives L := [[99, x, 0, 1], sin(x)].
\(L\) is a list. The second entry is \(f(x)\) and the first entry of \(L\) is a list which gives the coefficents of the ode. Notice they are ordered from lowest order to highest order of the ode. Since this is third order ode, then there are \(4\) entries. The first is the coefficient of \(y(x)\), the second is the coefficient of \(y'\) and the third is the coefficient of \(y''\) and the fourth entry is the coefficient of \(y'''\). Notice that coefficient of \(y''\) is zero since \(y''\) is missing from the ode.
This function DEtools:-convertAlg only works on linear ode’s. Therefore we need to
check if the ode is linear first before using. How to check for linear ode is given
above.
This function below finds the ode coefficients for any ode, linear or not. It takes as input the ode and returns back two items. The first is list of coefficients, and the second is the RHS (i.e. \(f(x)\) ).
So give an ode such as \( 3 y'' + 5 y' =\sin x\) it will return [0,5,3],sin(x) where the order of coefficents is \(y,y',y'',\dots \). If
terms is missing, zero is placed there.
#version Dec 3, 2024, 3PM. by Nasser M. Abbasi get_ode_coefficients_any_order:=proc(ode_in::`=`,func::function(name),$)::list,anything; #a,f #find coefficients. For example, for 5*y''''+y'''+y'+3*y=sin(x); it will return # [3,1,0,1,5],sin(x) #Where the list is in order y,y',y'',y''',y'''' local y:=op(0,func); local x:=op(1,func); local ode:=ode_in,i,j,tmp; local a::list; local c:=0; local f:=0; local ode_order::integer; local found_diff::truefalse; ode:= expand(lhs(ode))-expand(rhs(ode)); ode_order:= PDEtools:-difforder(ode,x); if ode_order<1 then error "ODE expected in get_ode_coefficients_any_order, found ", ode; fi; a:=[0$ode_order+1]; #for [y,y',y'',...,y(n)], that is why added 1 if ode::`+` then for i from 1 to nops(ode) do tmp := op(i,ode); found_diff:=false; for j from ode_order to 1 by -1 do if has(tmp, diff(y(x),x$j) ) then a[j+1] := a[j+1] + tmp/ diff(y(x),x$j); found_diff:=true; break; fi; od; if not found_diff then if has(tmp,y(x)) then a[1] := a[1] + tmp/y(x); else f := f+tmp; fi; fi; od; else found_diff:=false; for j from ode_order to 1 by -1 do if has(ode, diff(y(x),x$j) ) then a[j+1] := ode/ diff(y(x),x$j); found_diff:=true; break; fi; od; if not found_diff then if has(ode,y(x)) then a[1] := ode/y(x); else f :=ode; fi; fi; fi; return a,(-f); end proc:
Examples
ode:=sin(x)+y(x)*diff(y(x),x$5)+diff(y(x),x)-1/x+99*diff(y(x),x$3)=0; get_ode_coefficients_any_order(ode,y(x)) 1 [0, 1, 0, 99, 0, y(x)], -sin(x) + - x ode:=diff(y(x),x$3)=0: get_ode_coefficients_any_order(ode,y(x)) [0, 0, 0, 1], 0 ode:=diff(y(x),x$3)=sin(y(x)): get_ode_coefficients_any_order(ode,y(x)) [ sin(y(x)) ] [- ---------, 0, 0, 1], 0 [ y(x) ] ode:=sin(x)+1/x*diff(y(x),x$3)+y(x)=exp(x)+9; get_ode_coefficients_any_order(ode,y(x)) [ 1] [1, 0, 0, -], -sin(x) + exp(x) + 9 [ x]
I need to find the degree of the highest derivative in an ode.
For example, if the input is
Then the degree for highest derivative of \(y\) w.r.t. \(x\) is 4. For
It will be 2.
This function returns the order and degree of such term.
find_all_derivatives_of_specific_order:=proc(expr,y::symbol,x::symbol,N::posint)::set; local t1,t2; if not has(expr,diff(y(x),x$N)) then return {}; fi; t1 := identical(diff(y(x),x$N))^anything; t2 := identical(diff(y(x),x$N)); return indets[flat](expr,{t1,t2}); #MUST use flat end proc: get_order_and_degree_of_largest_derivative:=proc(expr,y::symbol,x::symbol)::integer,anything; local the_order,the_degree; local cand::set; local the_exponent; local item; if not has(expr,diff(y(x),x)) then return 0,0; fi; the_order := PDEtools:-difforder(expr,x); cand := find_all_derivatives_of_specific_order(expr,y,x,the_order); if nops(cand)=0 then the_degree := 0; else the_degree:=1; for item in cand do if type(item,`^`) then the_exponent := op(2,item); if type( the_exponent,symbol) or not type( the_exponent,numeric) then #assume larger the_degree :=the_exponent; else if type(the_degree,symbol) or not type( the_exponent,numeric) then next; else if op(2,item)>the_degree then the_degree := op(2,item); fi; fi; fi; else if type(the_degree,symbol) then next; else if the_degree=0 then the_degree := 1; fi; fi; fi; od; fi; return the_order,the_degree; end proc:
And now it can be called as follows
ode:=diff(y(x),x$2)+diff(y(x),x)+y(x)=0; get_order_and_degree_of_largest_derivative(lhs(ode),y,x) #2,1 ode:=diff(y(x),x$2)*diff(y(x),x$2)^4+diff(y(x),x)+y(x)=0; get_order_and_degree_of_largest_derivative(lhs(ode),y,x) #2,5 ode:=y(x)=0; get_order_and_degree_of_largest_derivative(lhs(ode),y,x) #0,0 ode:=diff(y(x),x)^n=0; get_order_and_degree_of_largest_derivative(lhs(ode),y,x) #1,n
Given
How to move all terms with derivative to LHS side and everything else to RHS?
ode:= a*diff(y(x),x)+b*diff(y(x),x$2)+x+cos(x)+y(x)+c*diff(y(x),x)^2=sin(x); ode:=lhs(ode)-rhs(ode); LHS,RHS:=selectremove(has,ode,'diff'); new_ode:=LHS=-RHS;
I had a need to find all derivatives of form diff(y( anything), anything) in an ode so to
check that \(y\) argument is not different among them.
For an example, given
the result should be
One issuse is how to check for diff and also check that the dependent variable is \(y\) so as not to
pick other dependent variables such as \(z\) in this example. This was done by converting diff to
D otherwise it will not work.
restart; expr:=a*diff(y(x),x$2)*diff(y(x),x$3)-sqrt(1+ diff(y(x),x$2)^2)+1/sin(diff(y(x),x$5))=diff(y(z),z)+exp(y(x)+diff(y(x),x))+diff(r(x),x); #this finds all derivatives list_of_diffs:=indets(expr,'satisfies'(s->op(0,s)='diff' and op([0,1],convert(s,D))=y)); #This finds all dependent variables list_of_diffs:=convert(list_of_diffs,list); list_of_indep_variables := map(x->PDEtools:-Library:-GetIndepVars(x)[-1],list_of_diffs); #this converts it to set. If the ODE is valid, then the list_of_indep_variables should #have one entry $x$ in it and nothing else. list_of_indep_variables := convert(ListTools:-Flatten(list_of_indep_variables),set); if nops(list_of_indep_variables)>1 then error( cat("Only one independent variable expected in differential form, found ", convert(list_of_indep_variables,string)) ); fi; if list_of_indep_variables[1]<>x then error( cat("Independent variable expected in differential form not same as independent variable in function given ", convert(list_of_indep_variables,string)) ); fi;
Another option instead of doing all the above is to do this
expr:=a*diff(y(x),x$2)*diff(y(x),x$3)-sqrt(1+ diff(y(x),x$2)^2)+1/sin(diff(y(x),x$5))=diff(y(z),z)+exp(y(x)+diff(y(x),x))+diff(r(x),x); try PDEtools:-Library:-GetDepVars([y(x)],expr); catch error "functions with name [y] having different dependency: [[y(x), y(z)]]" end try;
The function PDEtools:-Library:-GetDepVars([y(x)],expr) checks that only \(y(x)\)
dependency shows up. It throws an error otherwise. So if an error is thrown, then this means
\(y\) shows up with different independent variables.
Sometimes it is useful to invert an ode. i.e. make the independent variable the dependent variable, and the dependent variable the independent. For example, given
We want the ode to become
This can be done as follows
restart; ode:=1+ (x/y(x)-sin(y(x) ))*diff(y(x),x)=0; tr:={x=u(t),y(x)=t}; ode:=PDEtools:-dchange(tr,ode); ode:=eval(ode,[t=y,u=x]); ode:=simplify(ode);
In this case, we can get rid of the denomator, but this is a manual step for now.
ode:=numer(lhs(ode))=0;
The above can now be solved more easily for \(x(y)\) than solving the orignal ode for \(y(x)\).
For say Bessel ode of order zero:
eq:= x^2*diff(y(x),x$2)+x*diff(y(x),x)+x^2*y(x)=0; DEtools[indicialeq](eq,x,0,y(x)); #x^2 = 0
The third argument above is the singularity point of interest. So we have two roots, both zero. These are now used for finding the power series solution \(y(x)\) if needed.
Another example, is Bessel of order 1
eq:= x^2*diff(y(x),x$2)+x*diff(y(x),x)+(x^2-1)*y(x)=0; DEtools[indicialeq](eq,x,0,y(x)); #x^2-1 = 0
To write \(y'(x)=x\), one way is diff(y(x),x)=x and another is D(y)(x)=x. To write \(y''(x)=x\), one way is
diff(y(x),x$2)=x and another is (D@@2)(y)(x)=x.
To convert from one form to another use convert(eq,diff) or convert(eq,D)
to solve \(\frac {\partial u(x,t)}{\partial t}=k \frac {\partial ^2 u(x,t)}{\partial x^2}\) with homogeneous dirichlet boundary conditions \(u(0,t)=0,u(L,t)=0\) the commands are
restart; pde:=diff(u(x,t),t)=k*diff(u(x,t),x$2); bc:=u(0,t)=0,u(L,t)=0; sol:=pdsolve([pde,bc]) assuming 0<L:
Which gives
Which can be made more readable as follows
sol:=algsubs(_Z1=n,sol): sol:=algsubs(Pi*n/L=lambda(n),sol);
For homogeneous Neumann B.C., at \(x=0\), let \(\frac {\partial u}{\partial x}=0\) and at \(x=L\) let \(u(L,t)=0\), the solution it gives looks different than my hand solution
restart; pde:=diff(u(x,t),t)=k*diff(u(x,t),x$2); bc:=D[1](u)(0,t)=0,u(L,t)=0; pdsolve([pde,bc]) assuming 0<L;
It gives
I need to look more into the above and see if this comes out to be the same as my hand solution.
Another example, with initial conditions now given
restart; pde:=diff(u(x,t),t)=k*diff(u(x,t),x$2); bc:=D[1](u)(0,t)=0,u(L,t)=0; ic:=u(x,0)=f(x); sol:=pdsolve([pde,bc,ic],u(x,t)) assuming 0<L; sol1:=algsubs(_Z2=n,sol);
The result is
Another example
restart; pde:=diff(u(x,t),t)=k*diff(u(x,t),x$2); bc:=D[1](u)(0,t)=0,u(L,t)=0; ic:=u(x,0)=3*sin(Pi*x/L)-sin(3*Pi*x/L); sol:=pdsolve([pde,bc,ic],u(x,t)) assuming 0<L; sol1:=algsubs(_Z2=n,sol);
Another example
restart; pde:=diff(u(x,t),t)=k*diff(u(x,t),x$2); bc:=u(0,t)=0,u(L,t)=0; ic:=u(x,0)=3*sin(Pi*x/L)-sin(3*Pi*x/L); sol:=pdsolve([pde,bc,ic],u(x,t)) assuming 0<L;
The above answer seems wrong. There is not even a summation in it. It is different from my hand solution. Look more into it.
Add this
expr:=diff(y(x),x); Typesetting:-Settings(typesetprime=true, prime=x):
The above will display the expression as \(y'(x)\). To make it now show the \(x\) do
expr:=diff(y(x),x); Typesetting:-Settings(typesetprime=true, prime=x): Typesetting:-Suppress(y(x));
Now it will show the expression as just \(y'\). For all the above to work, make sure you have
Typesetting level set to Extended in the GUI.
This is done inside Tools->Options->Display menu.
To clear all the above Typesetting, do restart or do Typesetting:-Unsuppress(y(x))
The Maple syntax for seeting initial and boundary conditions is very confusing, as compared to Mathematica, which seems to me to be simpler. So I wrote this to remind me of the syntax each time.
For PDE, assuming dependent variable is \(u(x,t)\) then
| Conditions | Maple code |
| \(u(0,t)=0\) | u(0,t)=0 |
| \(\frac {\partial u}{\partial x}=0\) at \(x=0\) | D[1](u)(0,t)=0 |
| \(\frac {\partial ^2 u}{\partial x^2}=0\) at \(x=0\) | D[1,1](u)(0,t)=0 |
| \(\frac {\partial ^3 u}{\partial x^3}=0\) at \(x=0\) | D[1,1,1](u)(0,t)=0 |
| \(\frac {\partial u}{\partial t}=0\) at \(t=0\) | D[2](u)(x,0)=0 |
| \(\frac {\partial ^2 u}{\partial t^2}=0\) at \(t=0\) | D[2,2](u)(x,0)=0 |
| \(\frac {\partial ^3 u}{\partial t^3}=0\) at \(t=0\) | D[2,2,2](u)(x,0)=0 |
Notice the syntax for the last one above. It is (D[1]@@2)(u)(0,t)=0 and not
(D@@2)[1](u)(0,t)=0
For an ODE, assuming dependent variable is \(y(x)\) then the syntax is
| Conditions | Maple code |
| \(y(0)=0\) | y(0)=0 |
| \(\frac {dy}{dx}=0\) at \(x=0\) | D(y)(0)=0 |
| \(\frac {d^2 y}{d x^2}=0\) at \(x=0\) | (D@@2)(y)(0)=0 |
Use dsolve(ode,Lie)
To find symmetries, do
DEtools:-symgen(ode,y(x),HINT=[c__1+c__2*x+c__3*y,c__4+c__5*x+c__6*y])
or just
DEtools:-symgen(ode,y(x))
To debug it do
stopat(`ODEtools/symgen`); before calling dsolve or DEtools:-symgen
given an ode
do change of variable \(u(x)=y(x)^2\)
restart; ode:=y(x)=diff(y(x),x)^3*y(x)^2+2*x*diff(y(x),x); new_ode:=PDEtools:-dchange({y(x)=sqrt(u(x))},ode,{u});
given an ode
do change of variable \(t=\tau + \pi \)
restart; ode:=diff(y(t),t$2)+y(t)=2*t; PDEtools:-dchange({t=tau+Pi},ode,known={t},unknown={tau},params=Pi)
it is important to use params=Pi. Watch what happens if we do not do that
restart; ode:=diff(y(t),t$2)+y(t)=2*t; PDEtools:-dchange({t=tau+Pi},ode,known={t},unknown={tau});
Which is not what we want.
This verifies solution given in https://math.stackexchange.com/questions/3477732/can-t-see-that-an-ode-is-equivalent-to-a-bessel-equation Where a change of variables on
Was made using
To produce the ode
In Maple the above is done using
restart; ode := zeta^2*diff(eta(zeta),zeta$2) + zeta*diff(eta(zeta),zeta) - (zeta^2 + n^2)*eta(zeta) = 0; the_tr:={zeta=beta*x^gamma,eta(zeta)=y(x)/x^alpha}; PDEtools:-dchange(the_tr,ode,{y(x),x},'known'={eta(zeta)},'uknown'={y(x)},'params'={alpha,beta,gamma}); simplify(%); numer(lhs(%))=0; simplify(numer(lhs(%))/(x^(1-alpha)))=0; numer(lhs(%))=0; collect(%,[y(x),diff(y(x),x),diff(y(x),x$2)]);
Which gives
Here is another example. Here to make change of variables to polar coordinates by making \(x=r \cos \theta \) and \(y=r \sin \theta \) The ode is
In Maple
restart; ode := (y(x)-x*diff(y(x),x))/sqrt(1+ diff(y(x),x)^2) = x^2+y(x)^2; the_tr:={x=r(t)*cos(t),y(x)=r(t)*sin(t)}; PDEtools:-dchange(the_tr,ode,{r(t),t},'known'={y(x)},'uknown'={r(t)});
Which gives
Here is another example. Where we want to change \(R(r)\) to \(y(x)\) everywhere
restart; ode:=diff(R(r),r$2)+diff(R(r),r)+R(r)=0; the_tr:={r=x,R(r)=y(x)}; PDEtools:-dchange(the_tr,ode,{y(x),x},'known'={R(r),r},'uknown'={y(x),x});
The format of the transformation is old_independent_variable=new_independent_variable
and old_dependent_variable=new_dependent_variable
Make a phase plot of
The phase plot has \(x\) on the x axis and has \(y\) on the y axis. It shows the family of solutions for different initial conditions.
restart; ode:=diff(y(x), x) = sqrt(y(x)^2 - 1); DEtools:-DEplot(ode,y(x),x=-2..2,y=1..3)
To show specific solution curve that passes via some initial conditions such as \(y(0)=2\) then do
restart; ode:=diff(y(x), x) = sqrt(y(x)^2 - 1); DEtools:-DEplot(ode,y(x),x=-2..2,y=1..3,[y(0)=2])
Make a phase plot of
By plotting \(x(t)\) vs \(x'(t)\) without solving the ODE.
restart; alias(DS=DynamicSystems): ode := diff(x(t),t$2) +1/2*diff(x(t),t)+ x(t) = u(t); sys:=DS:-DiffEquation(ode,'outputvariable'=[x(t)],'inputvariable'=[u(t)]); sys0:=DS:-StateSpace(sys); eq1:=diff(x1(t),t)=sys0:-a[1,..].Vector([x1(t),x2(t)]); eq2:=diff(x2(t),t)=sys0:-a[2,..].Vector([x1(t),x2(t)]); DEtools:-DEplot([eq1,eq2],[x1(t),x2(t)],t=0..35,[[x1(0)=1,x2(0)=1]],x1=-2..2,x2=-2..2, numpoints=200, linecolor=black, axes=boxed);
Given ode where \(y\) is dependent variable and \(x\) is independent variable, how to normalize the ode so that all terms with \(y(x)\) are on left side and everything is on the RHS? THis makes it easier to see what the forcing function is. i.e. we want to write the ode as
Here is function to do this. This works for any ode or any expression.
export move_y_x_to_each_side:=proc(expr_in::`=`,y::symbol,x::symbol)::`=`; local expr:=expr_in; local stuff_with_y,stuff_not_with_y; expr := lhs(expr)-rhs(expr); if not expr::`+` then RETURN(expr_in); fi; expr := numer(normal(expr)); stuff_with_y,stuff_not_with_y := selectremove(has,expr,y); stuff_not_with_y := convert(stuff_not_with_y,diff); RETURN(stuff_with_y = -stuff_not_with_y); end proc:
Examples
ode:=diff(y(x),x$2)+x*y(x)+sin(x)-3=tan(x); expr:=move_y_x_to_each_side(ode,y,x); # diff(y(x), x, x) + x*y(x) = 3 - sin(x) + tan(x)
And
ode:=diff(y(x),x$2)+x*y(x)+sin(x)-3=tan(x)+y(x)+Pi; expr:=move_y_x_to_each_side(ode,y,x); # diff(y(x), x, x) + x*y(x) - y(x) = 3 - sin(x) + tan(x) + Pi
And
ode:=move_y_x_to_each_side(diff(x(t),t$2)/x(t) + diff(x(t),t)/x(t) - F(t)/x(t) =1,x,t) expr:=move_y_x_to_each_side(ode,x,t); # diff(x(t), t, t) + diff(x(t), t) - x(t) = F(t)
Sometimes I make solution for an ode with constants of integrations. These must be of form
_Cn where \(n\) is positive integer, such as _C1 or _C5. Now to integrate this solution
once more, need to come up with new constant that does not already show in the
solution.
To find current constants in an expression use this
sol:=x-_C1*x+_C3*x^2; indets(sol,And(symbol, suffixed(_C, nonnegint))); {_C1, _C3}
Therefore, one way to make a new constant of integration, is to find the largest numbered one and increase by one. Like this
restart; sol:=x-_C1*x+_C3*x^2; myconstants:=indets(sol,And(symbol, suffixed(_C, nonnegint))); map(X->String(X),myconstants); map(X->X[3..],%); map(X->:-parse(X),%); n:=max(%); new_constant:=_C||(n+1);
Sometimes it is usefull to know if ode is missing \(y(x)\) as this allows one to do the substitution \(y'=u\) and reduce the order of the ode.
This is a function which takes in an ode and returns true if it is missing \(y(x)\) or false otherwise.
#added 2/18/2024. checks if an ode is missing y export is_ode_missing_y:=proc(ode_in::`=`,y::symbol,x::symbol)::truefalse; local ode:=ode_in; local ode_order::posint := PDEtools:-difforder(ode_in); local N::posint; ode:=lhs(ode)-rhs(ode); ode:=convert(ode,D); for N from 1 to ode_order do ode:=eval(ode,[(D@@N)(y)(x)=Y||N]); od; RETURN(not has(ode,y(x))); end proc;
To use do
restart; ode:=diff(y(x),x$2)+diff(y(x),x$4)+x=0; is_ode_missing_y(ode,y,x); #true ode:=diff(y(x),x$2)+diff(y(x),x$4)+x=sin(x)+y(x); is_ode_missing_y(ode,y,x); #false
The general Riccati ode is
This can be converted to second order linear ode in \(u(x)\) using the transformation
Which results, after some simplifications in the ode
The following Maple code does the above
riccati_ode:= diff(y(x),x)= f__0(x)+f__1(x)*y(x)+f__2(x)*y(x)^2; Typesetting:-Unsuppress('all'); #always do this. Typesetting:-Settings(prime=x,'typesetprime'=true); #this says to use y'(x) instead of dy/dx Typesetting:-Suppress(u(x)); # this says to use y' and not y'(x) PDEtools:-dchange( {y(x)= -diff(u(x),x)/(f__2(x)*u(x))},riccati_ode,{u}); (lhs-rhs)(%); numer(normal(%)); collect(%,[diff(u(x),x),diff(u(x),x$2)]); map(X->X/f__2(x),%); new_U_ode:=-%=0;
Gives
new_U_ode := -(f__2(x)*f__1(x) + diff(f__2(x), x))*diff(u(x), x)/f__2(x) + f__0(x)*f__2(x)*u(x) + diff(u(x), x, x) = 0
Sometimes an ode can be solved much easier if we change the roles, and make the independent variable as the dependent variable.
For example, the ode \(y'(x)=\frac {1}{e^{4 y(x)} + 2 x}\) which is not easy to solve, becomes \(x'(y)=2 x(y) + e^{4 y}\) which is linear. After solving the simpler ode, we solve for \(y\) from the result which gives the solution to the original ode.
To do the conversion, use the command convert(ode, y_x). Below are examples
original_ode:= diff(y(x),x)=1/(exp(4*y(x))+2*x); sol_of_original := dsolve(ode); y(x) = ln(-c__1 - sqrt(c__1^2 + 2*x))/2, y(x) = ln(-c__1 + sqrt(c__1^2 + 2*x))/2 new_ode:= convert(original_ode, y_x); diff(x(y), y) = exp(4*y) + 2*x(y) sol_of_new := dsolve(new_ode); sol_of_new := eval(sol_of_new,x(y)=x); sol_of_new := y(x) = solve(sol_of_new,y); y(x) = (ln(-c__1 + sqrt(c__1^2 + 2*x))/2, ln(-c__1 - sqrt(c__1^2 + 2*x))/2)
We see that we get the same solution ofcourse, but solving ode after changing the roles is much easier now.
This is a standalone function that solves \(y'=a x^n + b y^2\). This is called the specific or reduced Riccati ode. The formulas used are from Eqworld ode0106 and Dr Dobrushkin web page
#solves y' = a*x^n + b*y^2, using standard formulas for n=-2 and n<>-2 special_riccati_dsolve := proc(n,a,b,func::function(name),$) local y:=op(0,func); local x:=op(1,func); local sol; local C1,C2; if not (n::'integer' or n::'fraction') then print("argument n must be integer or fraction"); RETURN(NULL); fi; if n=-2 then proc() local z,sol_z,lambda; try sol_z:=timelimit(30,[solve(b*z^2+z+a=0,z)]); if nops(sol_z)=0 then error "Unable to solve for root"; fi; catch: error "Unable to solve for root"; end try; z:=sol_z[1]; #pick one, any will work sol:=y(x)= z/x - x^(2*b*z)/( b*x/(2*b*z+1)*x^(2*b*z) + c__1); end proc(); else proc() local k,w; k:=1+n/2; C1:= c__1; C2:= c__2; try if evalb(a*b >0) then w:= sqrt(x) * ( C1*BesselJ(1/(2*k),1/k*sqrt(a*b)*x^k) + C2* BesselY(1/(2*k),1/k*sqrt(a*b)*x^k) ); else w:= sqrt(x) * ( C1*BesselI(1/(2*k),1/k*sqrt(-a*b)*x^k) + C2* BesselK(1/(2*k),1/k*sqrt(-a*b)*x^k) ); fi; catch: w:= sqrt(x) * ( C1*BesselJ(1/(2*k),1/k*sqrt(a*b)*x^k) + C2* BesselY(1/(2*k),1/k*sqrt(a*b)*x^k) ); end try; sol:= simplify(-1/b*diff(w,x)/w); sol:= y(x)=eval(sol,C2=1); end proc(); fi; RETURN(simplify(sol)); end proc:
Example usages
a:=1:n:=-2:b:=1: ode:=diff(y(x),x)=a*x^n+b*y(x)^2; sol:=special_riccati_dsolve(n,a,b,y(x)); sol := y(x) = ((-x^(sqrt(3)*I) - 3*c__1)*sqrt(3)*I + 3*x^(sqrt(3)*I) + 3*c__1)/(2*(x^(sqrt(3)*I)*sqrt(3)*I - 3*c__1)*x) a:=5:n:=-8/3:b:=3: ode:=diff(y(x),x)=a*x^n+b*y(x)^2; sol:=special_riccati_dsolve(n,a,b,y(x)); sol := y(x) = -sqrt(15)*((sqrt(15)*x - 15*c__1*x^(2/3) + c__1*x^(4/3)/3)*cos(3*sqrt(15)/x^(1/3)) + sin(3*sqrt(15)/x^(1/3))*(sqrt(15)*c__1*x + 15*x^(2/3) - x^(4/3)/3))/(((c__1*x^(1/3)*sqrt(15) + 45)*cos(3*sqrt(15)/x^(1/3)) + 45*(-sqrt(15)*x^(1/3)/45 + c__1)*sin(3*sqrt(15)/x^(1/3)))*x^2) a:=5:n:=5:b:=2: ode:=diff(y(x),x)=a*x^n+b*y(x)^2; sol:=special_riccati_dsolve(n,a,b,y(x)); sol := y(x) = -sqrt(10)*x^(5/2)*(BesselY(-6/7, (2*sqrt(10)*x^(7/2))/7) + BesselJ(-6/7, (2*sqrt(10)*x^(7/2))/7)*c__1)/(2*c__1*BesselJ(1/7, (2*sqrt(10)*x^(7/2))/7) + 2*BesselY(1/7, (2*sqrt(10)*x^(7/2))/7))
Ofcourse, dsolve can be used to solve these also.