There is as yet no general solution f for the equation

where f is a function dependent on t, n>1 is a positive integer, and a, b, and c are given constants unequal to 0 in F of characteristic zero. The design of the equation is by D. Bernoulli [2] with c=0, n =2 , as studied by him to predict the effect of smallpox vaccination. If c≠0,n=2 then (1.1) is named the Riccati equation with constant coefficients [20], and n=3 is a special case of the classical Abel equation discussed on p. 24 in [8]. The importance of this equation nowadays originates from the need in the physics of fluid dynamics to develop e.g. accurate weather prediction, neurological implants for deaf and blind people, or control fluids in nano capacities (for medicine administration through the brain-body blood barrier), etc. Fluid dynamics models have the format of partial differential equations (PDEs) of a high order, most often in the class of non-linear evolution Equations (NLEEs). Nature often has a complicated relationship between the phase speed of traveling waves (for instance tsunamis) and the parameters in the PDEs predicting them. The order of an NLEE becomes higher if more precision is required, e.g. dispersion effects are taken into account, but most often these model equations then become unsolvable. Kudryashov [12] designed the `Simplest Equation Method' (SEM) to reduce the order of intractable partial differential equations (PDEs) by assuming a series solution such that it solves a lower order Riccati equation.

An example of such use of a lower-order differential equation is in [9] for the solution of otherwise unsolvable double dispersion equations, i.e. the Sharma-Tasso-Olver (STO) equation. To glimpse other examples of the application of the work reported in this paper is the Benjamin-Ono PDEs in [23] which hit upon a Chini equation. Also, Chini’s method is still under study in applications [17] to Gross-Pitaevskii PDEs.

In general, the field of fluid dynamics progresses if the SEM becomes wider applicable by raising the order of the SEM. This paper gives an automated technique, such that high order solutions become feasible for the SEM.

The main result of this paper is to aid application fields with a computer-assisted method for symbolic differential solutions of (1.1), if a, b, and c are constants. To this end, we devise an algorithmic i.e. a constructive, method to find series solutions of any power n>1, with illustrative applications from [8], [18], [21]. There is as yet no general solution for this differential equation (1.1), if n>1 and a, b, c are general, i.e. non-constant.

Recall the formal derivative of power series over a field F with characteristic zero with

an element of the power series ring F[[t]], where t is an indeterminate over F. Then is the formal derivative f' of f defined as

The following lemma is known for a long time (see Formula 6.361 in [1], Formula 0.314 in [5], Theorem 1.6c in [7], [14] Ch. 17 (1st edn.), and Ch. 21 (2nd edn.). We bring it for reference.

Lemma 1. Let F be a field of characteristic zero and ∑ k=0∞Aktk a formal power series in F[[t]] with the indeterminate t. Then for each positive integer n,

where C0=(A0)n, and

Optimization of computation of power series expansions has been coined JCPMiller Algorithm in [6]. Later publications of the algorithm are in [7], [14], [22] though the algorithm has a long history. It was apparently known to Euler, according to Henrici [7], p. 65, but Wimp in [22] p. 2 refers to Lord Rayleigh [19] in 1910. Knopp [11], however, attributes it to Glaisher in 1875.

Now we prove the main theorem of this section:

Theorem 2. Let F be a field of characteristic zero and f=∑ k=0∞Aktk be a power series solution for the differential equation (1.1)

where n>1 is a positive integer, and a, b, c, and A₀ are given elements in F. If C0=(A0)n and m≥1,

then

and

Proof. Assume that a, b, c, and f(0)=A₀ are given elements in F and let f=∑ k=0∞Aktk be a power series solution for the differential equation (1.1). It is clear that f′=∑ k=0∞(k+1)Ak+1tk. On the other hand, by Lemma 1, a power series raised to the power n is calculated with the formula

where C0=(A0)n and

Finally, from the differential equation (1.1), we obtain

This completes the proof.

Let us recall the following Theorem 2.4.2 in [3]:

Theorem 3. (Existence and Uniqueness for First-Order Nonlinear Equations)

Let the functions f and ∂f/∂y be continuous in some rectangle n<t<β and γ<y<δ containing the point (t0,y0). Then, in some interval t0−h<t<t0+h contained in n<t<β there is a unique solution y=φ(t) of the initial value problem y′=f(t,y) and y(t0)=y0.

Corollary 4. Let a,b, and c be arbitrary real constants with a≠0 and let n>1 be integer. Then, there is a unique solution y=φ(t) satisfying equation (1.1)

with y(t0)=y0 on an open interval of real numbers including x₀.

Proof. Let us define f(t,y)=ayn+by+c. Then it is obvious that the functions f and ∂f/∂y are continuous and so by Theorem 3, the solution is unique.

The first recurrence in Theorem 2 expresses Cm,m=0,1,2,… in terms of a, b, c, and A_m. The A_m are external for this algorithm, i.e. 'global' in Maple's terminology. Rewriting this by a recursive program in the mathematical symbolic programming language Maple, we obtain the listing below of Miller's algorithm with the shortcut name C.

Maple's symbolic output of successive calls of Miller's algorithm C(m,4,A₀) for m=0,1,…,5 is

Remark 5. The proof of Theorem 2 also holds if the assumed solution is a formal Laurent series ∑ k=q∞Aktk, where q ≠ 0 is an integer number and may or may not be negative. A polynomial or power series over a field with a nonzero constant term, i.e. A₀ ≠ 0 is named `unit' in [7]. Note that the formal Laurent series over a field form a field again. The algorithm CLaurent is the Laurent extension of Miller's algorithm C at p. 55 in [7], made available in Remark 6 hereafter.

Remark 6. The recurrences (2.2) and (2.3) of Theorem 1 are such that each A_m can be expressed in terms of a, b, c, and A₀, for each m≥1. Both A and C do necessarily execute in this order: A0⇒C0⇒A1⇒C1⇒…⇒Am−1⇒Cm−1⇒Am. The full computation is given in Theorem 2. The first two steps are A0⇒C0 and subsequently C0⇒A1. Thereafter follow Am=(aCm−1+bAm−1)/m,for all m>1.

An algorithm - embedding previous tools - to find all symbolic solutions of the equation (2.4), whenever the constants a, b, c, and A₀ are given, is such that if these constants are numerically substituted, the algorithm solves correctly Bernoulli, Riccati, Abel as well as Chini's equations [8]. Because of this generality, we name our algorithm ABCR in alphabetic order of the names of the original scholars Abel, Bernoulli, Chini and Riccati. The body of the algorithm ABCR first sets the initial values A₀, A₁, C₁ and thereafter it alternates - as said - between updating C_k and A_k until the stop condition is satisfied.

This algorithm has to be changed slightly if non-units solutions are needed, i.e. the assumed power series f has A₀=0. Then the call to Miller's C algorithm needs to replaced by a call to the CLaurent algorithm, as follows.

Below we give applications by showing outputs of our ABCR algorithm.

Example 7. The solution of Bernoulli's equation c=0 in (2.4) is known in general

where C is a constant depending on a,b and a boundary condition y(0).

For the particular case y′(t)=(y(t))2+y(t),y(0)=1, let n=2,a=1,b=1,c=0 in (1.1), and collate the known solution

to the result below, by executing ABCR(1, 1, 0, 7, 2, 1): by Maple. This gives the coefficients of the series solution by seq(Aktk,k=0…7);. The output is the exact sequence of coefficients of expansion of the known solution.

Example 8. All test cases in the table below, taken from the literature, give corect results:

Example 9. Solution of an hitherto unknown Chini equation with n=9, a=1, b=1, c=1 in (1.1):

Example 10. Execution of the Laurent version of the algorithm

on a truncated non-unit series for a fourth order equation, for example with q=1, hence A₁ ≠ 0, and a1t+a2t2+…+a7t7.

The Claurent algorithm gives coefficients of Cm,m=0…7, as follows.

An advantage of ourmethod for equations with constant coefficients is the result of a series expansion without effort to analyse whether it has periodic solutions, or not. Our method grants solutions of all types. The series of solutions are obtained by the property discovered by Euler: a particular solution f₁ generates a series of solutions u via the substitution f = u + 1/f₁.

For the moltenslag problem in metallurgy - a special case of the celebrated NLEEs in fluid flow dynamics - no general solution for the slag temperature y(t) flowing out of the furnace was found thus far, see [13]. Our method solves this problem in full generality. The slag outflow temperature is fully known and predictable by our symbolic output. Our result is new to this metallurgy application field.

Also, by our method, we generate the needed expansion of the solution of Kudryashov's [12] SEM (Simplest Equation Method), to simplify otherwise intractable partial differential equations (PDEs) in the fluid dynamics field.

The technique is general, easy to implement via the two algorithms above, and yields exact expansion results.