The symmetry method is a powerful method to solve differential equations by finding invariant transformations (point symmetries) of the ODE. For first order ODEs, an integrating factor can be constructed from the symmetry and the ODE can then be solved. For first order ODEs, there is no closed form algorithm to find a symmetry in general, but there are algorithms for finding symmetries of a certain form. In Cheb-Terrab and Kolokolnikov [1], a very general algorithm is presented that searches for linear symmetries of the form
$[\xi = F(x),\eta=P(x)y+Q(x)]$.
In Cheb-Terrab and Roche [2], an algorithm is presented for finding several simpler symmetries of the form
(1) $[\xi=F(x)*G(y), \eta=0]$, $[\xi=0, \eta=F(x)*G(y)]$
(2) $[\xi=F(x)+G(y), \eta=0]$, $[\xi=0, \eta=F(x)+G(y)]$
(3) $[\xi=F(x), \eta=G(x)]$, $[\xi=F(y), \eta=G(y)]$
(4) $[\xi=F(x), \eta=G(y)]$, $[\xi=F(y), \eta=G(x)]$
(5) $[\xi=ax+by+c, \eta=fx+gy+h]$
We see that the linear symmetries presented in [1] are a generalization of the search for symmetries of the form (3) presented in [2].
These methods have been implemented in the function ode1solve and can be found in the file ode1_lie.mac. Additionally, several other methods have been implemented that search for symmetries of specific ODEs. These algorithms are very fast but also limited in applicability. The implemented methods are symmetry searches for ODEs of the type: linear, inverse-linear, separable, exact, Riccati, Bernoulli and Abel. These symmetries can also be found using the search for more general symmetries mentioned above, but these dedicated methods are usually much faster.
As a performance test, the database of Kamke's first order ODEs was used [3]. In the database are 367 ODEs of first order and first degree. A number of these ODEs (Cheb-Terrab and Kolokolnikov mention 20 of the first order - first degree) are unsolvable because they are too general (these are ODEs 47, 48, 50, 55, 56, 74, 79, 82, 202, 219, 234, 235, 237, 250, 253, 269, 331). 178 of the ODEs are of type separable, linear, inverse linear and Bernoulli and therefore considered straightforward to solve. They were not considered in the evaluation of the symmetry method in [2]. These equations (as well as all homogeneous and Ricatti and Abel odes with constant invariants) all have linear symmetries. Of the 246 ODEs of first order and first degree that are left, 133 were solved by the Cheb-Terrab and Roche method. With the method of Cheb-Terrab and Kolokolnikov, another 20 ODEs of first order and first degree with linear symmetries can be found.
Our implementation currently finds (nontrivial) symmetries for 298 of the 347 solvable ODEs of first order and first degree, which amounts to a success rate of 85%. Every nontrivial symmetry of a first order ODE directly leads to an integrating factor, and therefore to a solution of the ODE. Most ODEs for which no symmetries could be found are currently outside the scope of the method. They are Riccati or (nonconstant invariant) Abel ODEs with special function solutions like Bessel functions.
The file ode1_lie.mac has to be loaded before use. ode1_lie.mac depends on separable.mac, ode1_abel.mac and ode_extra.mac. A database of first order odes from the book of Kamke [3] is stored as a list with name kamke1 in the file kamke1_1.mac. We first load both these files (the kill(all) command will remove any user defined functions and variables. This is useful if you are running this notebook interactively).
In [ ]:
kill(all);batch("~/mathematics/maxima_files/kamke1_1.mac");batch("~/mathematics/maxima_files/ode1_lie.mac");
The code is very new and might produce a lot of output to help in the development. To reduce the amount of messages, we set the keyword DEBUGFLAG to 1:
In [4]:
DEBUGFLAG:1;
Out[4]:
Let's first solve a differential equation with a known solution. To solve the general first order linear ODE, simply call ode1solve,
In [5]:
ode: 'diff(y,x)=f(x)*y+g(x);
Out[5]:
In [6]:
ode1solve(ode,y,x);
Out[6]:
We recognize the general solution as found in many standard textbooks. When a solution has been found, the keyword method indicates which solver was used to find the symmetry, integrating method or solution.
In [7]:
method;
Out[7]:
The default behavior is to return only the solution, if it is found. You can control the output using the options returnSymmetries, returnIntegratingFactor and returnSolution.
In [8]:
ode1solve(ode,y,x,'returnSymmetries=true,'returnIntegratingFactor=true,'returnSolution=true);
Out[8]:
The output is a list $X,\mu,y$ containing the symmetry $X=[\xi,\eta]$, the integrating factor $\mu$ and the explicit solution of the ODE. The solver might return more than one solution. It is also possible to set the output behavior globally. For instance, if you want to return only the symmetries,
In [9]:
returnSymmetries:true;returnIntegratingFactor:false;returnSolution:false;
Out[9]:
Out[9]:
Out[9]:
In [10]:
X:ode1solve(ode,y,x);
Out[10]:
In [11]:
DEBUGFLAG:1;ode:kamke1[101];
Out[11]:
Out[11]:
The function isBernoulli(ode,y,x) will return the coefficients $b(x),n,a(x)$ of a Bernoulli ODE $y'+a(x)y=b(x)y^n$, or false if the ODE is not a Bernoulli ODE.
In [12]:
isBernoulli(ode,y,x);
Out[12]:
In [13]:
X:ode1solve(ode,y,x);
Out[13]:
We can test the symmetry by substituting it back into the determining equations. If $X=[\xi,\eta]$ is a symmetry of the ODE, the result is 0:
In [14]:
checkSymmetries(X,ode,y,x);
Out[14]:
An integrating factor for the ODE $y'=\frac{Q(x,y)}{P(x,y)}$can be obtained from the symmetry using $\mu=\frac{1}{\xi Q - \eta P}$. The correctness of the integrating factor can be checked in a similar way. Note that we previously set returnSymmetries to true and returnIntegratingFactor to false globally. We want to return the integrating factors for only this call to ode1solve:
In [15]:
mu:ode1solve(ode,y,x,'returnSymmetries=false,'returnIntegratingFactor=true,returnSolution=false);
Out[15]:
In [16]:
isIntegratingFactor(mu[1],ode,y,x);
Out[16]:
the solution can be checked by substituting it back into the ode:
In [17]:
sol:ode1solve(ode,y,x,'returnSymmetries=false,'returnIntegratingFactor=false,'returnSolution=true)[1];
Out[17]:
In [18]:
subst(sol,ode);
Out[18]:
In [19]:
ev(%,nouns);
Out[19]:
In [20]:
ratsimp(%);
Out[20]:
With the flag DEBUGFLAG you can control the amount of output printed on the screen. DEBUGFLAG can be set to a value between 0 and 5. The following rules are used for the level of detail of the output:
For most users, DEBUGFLAG is 1 or 2.
In [21]:
DEBUGFLAG:2;
Out[21]:
In [22]:
ode:kamke1[44];
Out[22]:
In [23]:
sol:ode1solve(ode,y,x,'returnSymmetries=true,'returnIntegratingFactor=true,'returnSolution=true);
Out[23]:
In this case, the solver returns two explicit solutions of the Bernoulli ODE. For many ODEs, the solution can not be written in explicit form. Sometimes, the implicit form is prefered over the explict form because the solution is much simpler in implicit form. With the keyword returnExplicit you can tell the solver to try to write the solution in explicit form or to return the implicit solution. An implicit solution of the first order ode is a first integral.
In [24]:
DEBUGFLAG:1;
sol:ode1solve(ode,y,x,'returnSymmetries=false,'returnIntegratingFactor=false,'returnSolution=true,'returnExplicit=false)[1];
Out[24]:
Out[24]:
In the book of Kamke, the solution is showns semi-explicitly, with $1/y^2$ on the left-hand side:
In [25]:
ratexpand(solve(subst(y^2=1/Y,%),Y));
Out[25]:
In [26]:
ode:'diff(y,x)=f(x)*g(y);
Out[26]:
In [27]:
ode1solve(ode,y,x,'useMethod="separable",'returnSymmetries=false,'returnIntegratingFactor=false,'returnSolution=true)[1];
Out[27]:
The solution is returned in implicit form, because the solver failed to find the explicit solution. We get a warning because returnExplicit is still true. useMethod can also be a list of methods and the methods will be tried in the order as they appear in the list.
Using boundary conditions works in the same way as for the built-in ode solvers. We use ic1(solution,x=x0,y=y0) to impose the boundary condition y(x0)=y0 on the solution of the ode. Note that ode1solve takes into account the global values integration_constant and integration_constant_counter. The default integration constant for first order ODEs is %c as is the case with the built-in ode solver.
In [28]:
ode:'diff(y,x)=x;
Out[28]:
In [29]:
sol:ode1solve(ode,y,x,'returnSymmetries=false,'returnIntegratingFactor=false,'returnSolution=true)[1];
Out[29]:
In [30]:
FIX_INTEGRATION_CONSTANT:false;
Out[30]:
In [31]:
sol:ode1solve(ode,y,x,'returnSymmetries=false,'returnIntegratingFactor=false,'returnSolution=true)[1];
Out[31]:
In [32]:
sol:ode1solve(ode,y,x,'returnSymmetries=false,'returnIntegratingFactor=false,'returnSolution=true)[1];
Out[32]:
In [33]:
ic1(sol,x=x0,y=y0);
Out[33]:
In [34]:
ic1(sol,x=0,y=1);
Out[34]:
[1] E.S. Cheb-Terrab and T. Kolokolnikov, First-order ordinary differential equations, symmetries and linear transformations, Euro. J. of Applied Mathematics 14 (2003)
[2] E.S. Cheb-Terrab and A.D. Roche, Symmetries and first order ODE patterns, Computer Physics Communications 113 (1998)
[3] E. Kamke, Differentialgleichungen, L$\ddot o$sungsmethoden und L$\ddot o$sungen, Leipzig (1959)
[4] E.S. Cheb-Terrab, L.G.S. Duarte and L.A.C.P. da Mota, Computer algebra solving of first order ODEs using symmetry methods, Computer Physics Communications 101 (1997)
[5] F. Schwarz, Symmetry analysis of Abel's equation, Studies in applied mathematics 100 (1998)
[6] F. Schwarz, Algorithmic solution of Abel's equation, Computing 61 (1998)