Diﬀerential equations

3 Diﬀerential equations

3.1 How to check if ode is of certain type?
3.2 How to force dsolve to use speciﬁc method for solving?
3.3 How to ﬁnd a particular solution to ODE?
3.4 How to ﬁnd basis solutions for homogeneous ode?
3.5 How to solve a diﬀerential equation with initial conditions?
3.6 How to verify that the ODE solution given is correct?
3.7 How to know the type of ODE?
3.8 What packages to load for diﬀerential equations?
3.9 How to plot solution of diﬀerential equations?
3.10 On High precision. Using taylor to solve ODE
3.11 Obtain ODE in canonical coordinates (Lie symmetry)
3.12 How to use hint with symgen?
3.13 How to parse a single ode?
3.14 How to check if single ODE is valid?
3.15 How to check if an ode has $y'$ in it?
3.16 How to check if an ode is linear ode?
3.17 How to ﬁnd the order of an ode?
3.18 How to ﬁnd the coeﬃcients of a linear ode?
3.19 ﬁnd order and degree of highest derivative
3.20 How to move all derivatives to one side in an equation?
3.21 How to obtains list of all derivatives in expression?
3.22 How to invert roles of dependent variable and independent variable in an ode?
3.23 How to ﬁnd the indicial equation for an ODE?
3.24 How to write derivative
3.25 How to solve heat PDE in 1D in Maple 2017?
3.26 How to make Maple display diﬀ(y(x),x) as y’(x) or as y’ ?
3.27 How to set boundary conditions for dsolve or pdsolve?
3.28 How force dsolve to use Lie?
3.29 How to do change of variables on the dependent variable for an ODE?
3.30 How to do change of variable on the independent variable for an ODE?
3.31 ODE change of variable on both dependent and independent variable?
3.32 How to make phase plot of ﬁrst order ODE?
3.33 How to make phase plot of second order ODE?
3.34 How to move all terms with y to one side?
3.35 How to ﬁnd all constants of integrations in an expression?
3.36 How to ﬁnd if an ode is missing the dependent variable?

3.1 How to check if ode is of certain type?

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 speciﬁc 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]

3.2 How to force dsolve to use speciﬁc method for solving?

To ﬁnd 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 ﬁrst 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 ﬁrst 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 speciﬁc 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 ﬁnd all methods, and for higher order ode, we ﬁrst 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 ﬁnd 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]

3.3 How to ﬁnd a particular solution to ODE?

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`);

3.4 How to ﬁnd basis solutions for homogeneous ode?

Use the output=basis option

ode:=diff(y(x),x$2)-x*diff(y(x),x)-x*y(x)=0; 
dsolve(ode,output=basis);

3.5 How to solve a diﬀerential equation with initial conditions?

To solve

\[ y''-3y'+2y=10 e^{5 x} \]

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));

3.6 How to verify that the ODE solution given is correct?

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

3.7 How to know the type of ODE?

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');

3.8 What packages to load for diﬀerential equations?

Use with(DEtools);

3.9 How to plot solution of diﬀerential equations?

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);

3.10 On High precision. Using taylor to solve ODE

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.

3.11 Obtain ODE in canonical coordinates (Lie symmetry)

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 ﬁrst 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 ﬁrst 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 ﬁrst one obtained by Maple. This makes it possible to experiment with diﬀerent inﬁnitesimals. Note that diﬀerent inﬁnitesimals will/can result in diﬀerent canonical ODE, but the ﬁnal 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)

3.12 How to use hint with symgen?

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]

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.

3.13 How to parse a single ode?

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));

3.14 How to check if single ODE is valid?

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"

3.15 How to check if an ode has $y'$ in it?

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 ﬁrst 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 ﬁrst order diﬀ 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

3.16 How to check if an ode is linear ode?

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.

3.17 How to ﬁnd the order of an ode?

ode:=diff(y(x),x$3)+ x*diff(y(x),x)+sin(x)=0; 
PDEtools:-difforder(ode)

Gives $3$.

3.18 How to ﬁnd the coeﬃcients of a linear ode?

Given a linear ode of the form $A y''+B y' + C y= f(x)$ how to ﬁnd $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 ﬁrst entry of $L$ is a list which gives the coeﬃcents 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 ﬁrst is the coeﬃcient of $y(x)$, the second is the coeﬃcient of $y'$ and the third is the coeﬃcient of $y''$ and the fourth entry is the coeﬃcient of $y'''$. Notice that coeﬃcient 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 ﬁrst before using. How to check for linear ode is given above.

3.19 ﬁnd order and degree of highest derivative

I need to ﬁnd the degree of the highest derivative in an ode.

For example, if the input is

\[ x \left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right )^{5}+x \left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right )^{2}+\left (\frac {d}{d x}y \left (x \right )\right ) y \left (x \right )+\left (\frac {d^{3}}{d x^{3}}y \left (x \right )\right )^{4} \sin \left (x \right )+5 y \left (x \right )+\sin \left (x \right )+\frac {d^{6}}{d x^{6}}r \left (x \right )+\left (\frac {d^{3}}{d x^{3}}y \left (x \right )\right )^{4} \cos \left (x \right ) = 0 \]

Then the degree for highest derivative of $y$ w.r.t. $x$ is 4. For

\[ y''(x)^2 + y'(x)+y(x)=0 \]

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

3.20 How to move all derivatives to one side in an equation?

Given

\[ a \left (\frac {d}{d x}y \left (x \right )\right )+b \left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right )+x +\cos \left (x \right )+y \left (x \right )+c \left (\frac {d}{d x}y \left (x \right )\right )^{2} = \sin \left (x \right ) \]

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;

\[ a \left (\frac {d}{d x}y \left (x \right )\right )+b \left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right )+c \left (\frac {d}{d x}y \left (x \right )\right )^{2} = -x -\cos \left (x \right )-y \left (x \right )+\sin \left (x \right ) \]

3.21 How to obtains list of all derivatives in expression?

I had a need to ﬁnd all derivatives of form diff(y( anything), anything) in an ode so to check that $y$ argument is not diﬀerent among them.

For an example, given

\[ a \left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right ) \left (\frac {d^{3}}{d x^{3}}y \left (x \right )\right )-\sqrt {1+\left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right )^{2}}+\frac {1}{\sin \left (\frac {d^{5}}{d x^{5}}y \left (x \right )\right )} = \frac {d}{d z}y \left (z \right )+{\mathrm e}^{y \left (x \right )+\frac {d}{d x}y \left (x \right )}+\frac {d}{d x}r \left (x \right ) \]

the result should be

\[ \left [\begin {array}{c} \frac {d^{5}}{d x^{5}}y \left (x \right ) \\ \frac {d^{4}}{d x^{4}}y \left (x \right ) \\ \frac {d^{3}}{d x^{3}}y \left (x \right ) \\ \frac {d^{2}}{d x^{2}}y \left (x \right ) \\ \frac {d}{d x}y \left (x \right ) \\ \frac {d}{d z}y \left (z \right ) \end {array}\right ] \]

One issuse is how to check for diﬀ 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 diﬀerent independent variables.

3.22 How to invert roles of dependent variable and independent variable in an ode?

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

\[ 1+\left (\frac {x}{y \left (x \right )}-\sin \left (y \left (x \right )\right )\right ) \left (\frac {d}{d x}y \left (x \right )\right ) = 0 \]

We want the ode to become

\[ -\sin \left (y \right ) y +y \left (\frac {d}{d y}x \left (y \right )\right )+x \left (y \right ) = 0 \]

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);

\[ \frac {-\sin \left (y \right ) y +y \left (\frac {d}{d y}x \left (y \right )\right )+x \left (y \right )}{y \left (\frac {d}{d y}x \left (y \right )\right )} = 0 \]

In this case, we can get rid of the denomator, but this is a manual step for now.

ode:=numer(lhs(ode))=0;

\[ -\sin \left (y \right ) y +y \left (\frac {d}{d y}x \left (y \right )\right )+x \left (y \right ) = 0 \]

The above can now be solved more easily for $x(y)$ than solving the orignal ode for $y(x)$.

3.23 How to ﬁnd the indicial equation for an ODE?

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 ﬁnding 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

3.24 How to write derivative

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)

3.25 How to solve heat PDE in 1D in Maple 2017?

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

\[ u \left ( x,t \right ) =\sum _{{\it \_Z1}=1}^{\infty }{\it \_C1} \left ( {\it \_Z1} \right ) \sin \left ( {\frac {\pi \,{\it \_Z1}\,x}{L}} \right ) {{\rm e}^{-{\frac {k{\pi }^{2}{{\it \_Z1}}^{2}t}{{L}^{2}}}}} \]

Which can be made more readable as follows

sol:=algsubs(_Z1=n,sol): 
sol:=algsubs(Pi*n/L=lambda(n),sol);

\[ u \left ( x,t \right ) =\sum _{n=1}^{\infty }{\it \_C1} \left ( n \right ) \sin \left ( x\lambda \left ( n \right ) \right ) {{\rm e}^{-kt \left ( \lambda \left ( n \right ) \right ) ^{2}}} \]

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 diﬀerent 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

\[ u \left ( x,t \right ) ={\it \_C3}\,{\it \_C2}\, \left ( {{\rm e}^{1/4\,{\frac {2\,i\pi \,xL-k{\pi }^{2}t}{{L}^{2}}}}}+{{\rm e}^{-1/4\,{\frac {\pi \, \left ( 2\,ixL+k\pi \,t \right ) }{{L}^{2}}}}} \right ) \]

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

\[ u \left ( x,t \right ) =\sum _{n=1}^{\infty } \left ( 2\,{\frac {1}{L}{{\rm e}^{-1/4\,{\frac {k{\pi }^{2}t \left ( 1+2\,n \right ) ^{2}}{{L}^{2}}}}}\cos \left ( 1/2\,{\frac {\pi \,x \left ( 1+2\,n \right ) }{L}} \right ) \int _{0}^{L}f \left ( x \right ) \cos \left ( 1/2\,{\frac {\pi \,x \left ( 1+2\,n \right ) }{L}} \right ) \,{\rm d}x} \right ) \]

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);

\[ u \left ( x,t \right ) =\sum _{n=1}^{\infty }768\,{\frac {1}{\pi \, \left ( 16\,{n}^{4}+32\,{n}^{3}-136\,{n}^{2}-152\,n+105 \right ) }{{\rm e}^{-1/4\,{\frac {k{\pi }^{2}t \left ( 1+2\,n \right ) ^{2}}{{L}^{2}}}}}\cos \left ( 1/2\,{\frac {\pi \,x \left ( 1+2\,n \right ) }{L}} \right ) } \]

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;

\[ u \left ( x,t \right ) =\sin \left ( {\frac {\pi \,x}{L}} \right ) {{\rm e}^{-9\,{\frac {{\pi }^{2}kt}{{L}^{2}}}}} \left ( -2\,\cos \left ( 2\,{\frac {\pi \,x}{L}} \right ) +3\,{{\rm e}^{8\,{\frac {{\pi }^{2}kt}{{L}^{2}}}}}-1 \right ) \]

The above answer seems wrong. There is not even a summation in it. It is diﬀerent from my hand solution. Look more into it.

3.26 How to make Maple display diﬀ(y(x),x) as y’(x) or as y’ ?

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))

3.27 How to set boundary conditions for dsolve or pdsolve?

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`

3.28 How force dsolve to use Lie?

Use dsolve(ode,Lie)

To ﬁnd 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

3.29 How to do change of variables on the dependent variable for an ODE?

given an ode

\[ y \left (x \right ) = \left (\frac {d}{d x}y \left (x \right )\right )^{3} y \left (x \right )^{2}+2 x \left (\frac {d}{d x}y \left (x \right )\right ) \]

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});

\[ \sqrt {u \left (x \right )} = \frac {\left (\frac {d}{d x}u \left (x \right )\right )^{3}}{8 \sqrt {u \left (x \right )}}+\frac {x \left (\frac {d}{d x}u \left (x \right )\right )}{\sqrt {u \left (x \right )}} \]

3.30 How to do change of variable on the independent variable for an ODE?

given an ode

\[ \frac {d^{2}}{d t^{2}}y \left (t \right )+y \left (t \right ) = 2 t \]

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)

\[ \frac {d^{2}}{d \tau ^{2}}y \left (\tau \right )+y \left (\tau \right ) = 2 \tau +2 \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});

\[ \frac {d^{2}}{d \tau ^{2}}y \left (\tau , \pi \right )+y \left (\tau , \pi \right ) = 2 \tau +2 \pi \]

Which is not what we want.

3.31 ODE change of variable on both dependent and independent variable?

This veriﬁes 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

\[ \xi ^2 \frac {d^2\eta }{d\xi ^2} + \xi \frac {d\eta }{d\xi } - (\xi ^2+n^2)\eta =0 \]

Was made using

\[ \eta =\frac {y}{x^\alpha }, \quad \xi =\beta x^\gamma , \]

To produce the ode

\[ \frac {d^2y}{dx^2} - \frac {(2\alpha -1)}{x}\frac {dy}{dx} - (\beta ^2 \gamma ^2 x^{2\gamma -2}+\frac {n^2\gamma ^2-\alpha ^2}{x^2})y=0. \]

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

\[ \left (-\gamma ^{2} x^{-1+2 \gamma } \beta ^{2} x -\gamma ^{2} n^{2}+\alpha ^{2}\right ) y \left (x \right )+\left (-2 \alpha x +x \right ) \left (\frac {d}{d x}y \left (x \right )\right )+x^{2} \left (\frac {d^{2}}{d x^{2}}y \left (x \right )\right ) = 0 \]

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

\[ \frac {y-x y'}{\sqrt {1+(y')^2}}=x^2+y^2 \]

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

\[ \frac {r \left (t \right ) \sin \left (t \right )-\frac {r \left (t \right ) \cos \left (t \right ) \left (\left (\frac {d}{d t}r \left (t \right )\right ) \sin \left (t \right )+r \left (t \right ) \cos \left (t \right )\right )}{\left (\frac {d}{d t}r \left (t \right )\right ) \cos \left (t \right )-r \left (t \right ) \sin \left (t \right )}}{\sqrt {1+\frac {\left (\left (\frac {d}{d t}r \left (t \right )\right ) \sin \left (t \right )+r \left (t \right ) \cos \left (t \right )\right )^{2}}{\left (\left (\frac {d}{d t}r \left (t \right )\right ) \cos \left (t \right )-r \left (t \right ) \sin \left (t \right )\right )^{2}}}} = r \left (t \right )^{2} \left (\cos ^{2}\left (t \right )\right )+r \left (t \right )^{2} \left (\sin ^{2}\left (t \right )\right ) \]

Here is another example. Where we want to change $R(r)$ to $y(x)$ everywhere

\[ \frac {d^{2}}{d r^{2}}R \left (r \right )+\frac {d}{d r}R \left (r \right )+R \left (r \right ) = 0 \]

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});

\[ \frac {d^{2}}{d x^{2}}y \left (x \right )+\frac {d}{d x}y \left (x \right )+y \! \left (x \right ) = 0 \]

The format of the transformation is old_independent_variable=new_independent_variable and old_dependent_variable=new_dependent_variable

3.32 How to make phase plot of ﬁrst order ODE?

Make a phase plot of

\[ y'(x) = \sqrt {y(x)^{2}-1} \]

The phase plot has $x$ on the x axis and has $y$ on the y axis. It shows the family of solutions for diﬀerent 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 speciﬁc 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])

3.33 How to make phase plot of second order ODE?

Make a phase plot of

\[ \frac {d^{2}}{d t^{2}}x \left (t \right )+\frac {\frac {d}{d t}x \left (t \right )}{2}+x \left (t \right ) = u \left (t \right ) \]

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);

3.34 How to move all terms with y to one side?

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

\[ a y''+b y'+c y= f(x) \]

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)

3.35 How to ﬁnd all constants of integrations in an expression?

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 ﬁnd 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 ﬁnd 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);

3.36 How to ﬁnd if an ode is missing the dependent variable?

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

[next] [prev] [prev-tail] [front] [up]


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`


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`

3 Diﬀerential equations

3.1 How to check if ode is of certain type?

3.2 How to force dsolve to use speciﬁc method for solving?

3.3 How to ﬁnd a particular solution to ODE?

3.4 How to ﬁnd basis solutions for homogeneous ode?

3.5 How to solve a diﬀerential equation with initial conditions?

3.6 How to verify that the ODE solution given is correct?

3.7 How to know the type of ODE?

3.8 What packages to load for diﬀerential equations?

3.9 How to plot solution of diﬀerential equations?

3.10 On High precision. Using taylor to solve ODE

3.11 Obtain ODE in canonical coordinates (Lie symmetry)

3.12 How to use hint with symgen?

3.13 How to parse a single ode?

3.14 How to check if single ODE is valid?

3.15 How to check if an ode has \(y'\) in it?

3.16 How to check if an ode is linear ode?

3.17 How to ﬁnd the order of an ode?

3.18 How to ﬁnd the coeﬃcients of a linear ode?

3.19 ﬁnd order and degree of highest derivative

3.20 How to move all derivatives to one side in an equation?

3.21 How to obtains list of all derivatives in expression?

3.22 How to invert roles of dependent variable and independent variable in an ode?

3.23 How to ﬁnd the indicial equation for an ODE?

3.24 How to write derivative

3.25 How to solve heat PDE in 1D in Maple 2017?

3.26 How to make Maple display diﬀ(y(x),x) as y’(x) or as y’ ?

3.27 How to set boundary conditions for dsolve or pdsolve?

3.28 How force dsolve to use Lie?

3.29 How to do change of variables on the dependent variable for an ODE?

3.30 How to do change of variable on the independent variable for an ODE?

3.31 ODE change of variable on both dependent and independent variable?

3.32 How to make phase plot of ﬁrst order ODE?

3.33 How to make phase plot of second order ODE?

3.34 How to move all terms with y to one side?

3.35 How to ﬁnd all constants of integrations in an expression?

3.36 How to ﬁnd if an ode is missing the dependent variable?