Most physical and chemical processes encountered in chemical and process engineering often lead to differential equations which are as varied as the complex mass, momentum and energy transport processes they represent. Thus finding/selecting the appropriate one to use in any given situation is not an easy task. Following the realization that the solution of any differential equation can be generally expressed as a polynomial (truncated power series) which can be regressed by the least square method, and the coefficients of the regressed model linked to that of the binomial formula led to a generalized computational procedure for solving a wide-range of differential equations in chemical and process engineering (both initial and boundary value problems). The method was found to be computationally efficient and inexpensive as it is fast converging and free from rounding off errors and overshoot.

Keywords: Physical and chemical processes; Differential equations, initial and boundary value problems; least square regression analysis; binomial formula; computational procedure

Differential equations-used in most science and engineering analysis- are particularly important in modeling chemical engineering systems. As the physical and chemical laws governing these processes (heat, mass, and momentum transfer) are complex, as well as the chemical reactions, reaction heat, adsorption, desorption, phase transition, multiphase flow, etc,associated with them, so also are the differential models arising from them and include:, homogeneous/non-homogeneous, linear/non-linear , constant/variable coefficients, 1^st, 2^nd or higher orders, and systems of simultaneous differential equationsboth for Ordinary and partial differential equations. Numerous methods have been developed for solving them. They include:

Exact methods such as Method of undetermined coefficients, Integrating factor, Method of variation of parameters/separation of variables[1,2].

Approximate (but convergent) methods such as Successive Approximations, Perturbation Theory, Multiple scale Analysis, power series solutions, Generalized Fourier series/orthogonal functions [3].

Numerical methods such as: Euler methods; Trapezoidal rule; Runge–Kutta methods; Finite difference methods ;Finite element methods, gradient discretisationmethods[4-18].

Numerical techniques are particularly important since most realistic deferential equations in chemical and process engineering do not have exact analytic solutions, therefore, approximate/numerical techniques, such as Runge-Kutta, Euler, Newton gradient methods, finite difference/element) are used extensively.These methods have proven rather successful in dealing with both linear and nonlinear problems as well as initial boundary problems(IVP) and boundary value problems(BVP). Despite these obvious advantages, they have some demerits. Since they involve discretization, they have rounding off errors and are computationally expensive, and in some cases will not converge to the true solution.

particular generalized procedure which is used to grind most differential equations is the Power Series Method(PSM) (https://www.researchgate.net/publication/293652327). This method has proven rather successfulin dealing with both linear as well as nonlinear problems, as it yields analytical solutions andoffers certain advantages over standard numerical methods. It is free from rounding off errorssince it does not involve discretization, and is computationally inexpensive. In spite of the advantages of PSM, it has some drawbacks, the major one being that it cannot handle boundary value problems, as it can only handle initial value problems.Thus finding/selecting the appropriate one to use in any given situation is not an easy task.

In this article, we present a more generalized computational approach, which is amodification of PSMusing the Binomial formula/least square approximation, the goal being to use the same platform/ procedure tosimulate all differential equations in chemical and process engineering(both initial and boundary value problems),and to overcome the shortcoming of other methods.

Usually we solve differential equations to obtain the value of the dependent variable in terms of independent variable/s, which are finally represented in graphs for better visualization and analysis. These solution curves are as varied as the differential equations they are representing. Thus the first step in unification/generalization is to find a general expression that can fit all data/graphs as accurately as possible. We know from statistics that any set of data (or expression) can be regressed by the general least square method [19,20], and generally, the higher the order of the regression, the more accurate the solution becomes.Thusto fit an appropriate equation to the data, we use regression analysis given generally as:

Y = a0 + a1X + a2X² + a3X³ + … + anXⁿ(1.1)

The coefficients , ai , are found by least square method, applied to the X-Y data,which results in the matrix:

Similarly the solution of any differential equation can be expressed and regressed as a polynomial of the eq.1 form.Notice that polynomials are simply finite or truncated power series. That is a polynomial is a power series where the coefficients, ai beyond a certain point are all zero. Thus the coefficients of the polynomial can be found by manipulating power series by differentiation and theuse the Taylor series expansion. These manipulations will lead to a recurrent relation- recursion equation- given as:

The above equation is called a recurrence relation for the coefficients of the power series. For instance the first three equations of the relation are:

Hence, once the first two coefficients, a0 and a1, are known(which are actually the initial boundary conditions) then all other coefficients are determined by the recurrence relation. Thus PSM is limited to only initial boundary value (IVP) problems as it cannot handle boundary value problems (BVP) because of its recursive nature. In a nutshell, the method admits a polynomial expression as the solution of the differential equation, which when differentiated and substituted into the differential equation gives a recurrence relation-recursive formula, which by knowing the initial boundary conditions, gives the coefficients of the polynomial expression, hence the solution of the differential equation. This approach is adopted, with some modifications, in solving differential equations(both initial and boundary value problems). All that is needed is to find the right order and the right parameters(coefficients) of the polynomial expression(representing the differential equation); by modeling/numericalanalysis as this work shows.

Let us start with a general second order linear ODE with both constant and variable coefficients as follows:

f_i(t) can be any function of t namely: linear, nonlinear, exponential, trigonometric,Logarithmic, or their combinations. Whatever the form or forms of f(t), the general procedure is to remodel and express in terms of polynomial expressions within the solution interval, using the least square regression analysis already explained. Thus:

therefore t - t₀ = t- t₁ = t- t₂ = t- t₃ = … = t-t_n = 0 (12) This implies that: t –t_i = 0_i= 0,1,2,3…,n(13) Thus (t - t_i)ⁿ = 0t₀≤t_i≤t_n (14)

Eqs. 10 and 14 are equivalent. This realization/discovery that eqs.10 and 14 are equivalent is the key to new computational method for generalization of the method of simulating differential equations (both IVP and BVP). Thus the parameters or coefficients of eq,14 are linked to those of eq.10using binomial formula:

Thus eq.17 and eq.10 are equivalent, and term by term comparison result in the following equations that linked their coefficients:

Thus there are n-equations (eqs. 18.1 to 18.n) with n+1 unknown parameters, a₀,a₁,a₂,a₃,a₄,a₅,…,a_n thereby requiring additional equation/s, which are found from the boundary conditions of the problem. Usually the boundary conditions can be one( for IVP of first order) or two ( for BVP of first order and IVP/BVP of second order). Let us consider the two cases, and clearly visualize the procedures by considering n=5. Thus for n=5, eq.18 set reduces to:

Thus we have five equations (eq.19.1 to 19.5) with six unknowns (a0 to a6). The remaining equation/s will come from the boundary conditions as follows:

Case 1: One point boundary value problem. for 1st order IVP, whereY(τ) is given, thereby giving one boundary equation, eq.4 in terms of the boundary condition provides the remaining equation(for n=5) as:

Y(τ) = a0 + a1τ + a2τ² + a3τ³+ a4τ⁴ + a5τ⁵ (19.6)

Rearrangingeqs. 19.1 to 19.6 result in the following sample simulation matrix for one point IVP ODE’s:

Case 2: Two points boundary value problem.

Similarly for BVP of first order ODE or both for IVP and BVP of second order ODE,s where Y(τ) and Y’(τ)or Y(𝞽1)and Y’(𝞽2) or Y(τ1) and Y(τ2) pairsare given, thereby given two boundary equations in each case as follows:

This implies we now have five-main equations(eqs.19.1 to 19.5) and two boundary equations (eqs. 21.1 and 21.2 or 22.1 and 22.2 or 23.1 and 23.2) with only six parameters (a0,a1,a2,a3,a4 and a5). Thus we need to reduce the main equations from n=5 to n-1=4 to accommodate the two equations from the boundary conditions. Thus eqs.19.1 to 19.4 are divided by eq.19.5 to give the required four equations as follows:

Rearranging and combining eqs.24.1 to 24.4 with the boundary value equations result in the following sample simulation matrix in each case as follows:

A closer observation on the sample simulation matrices show certain similar trends that can lead to generalization and extension to higher order and matrix size. Thus a spreadsheet program has been developed by the authors for handling higher order ODE’s and matrix size.

The analyses here follow the same pattern as the linear case explained in the previous section. Thus let us illustrate by converting the general second order linear ODE of eq3 into second order non linear form by making the y-coefficientvariable in the dependent variable as highlighted in eq.26:

(k0+k1t + k2t2)d2Y/ dt2 + (k3+k4t) dY/dt + (k5 + k6Y) Y + f(t) + k7 = 0 (26)

Eq.26 is a second order nonlinear ODE. Thus by following similar regression procedure as explained in the linear cases, the resulting equation(for the sample case of n=5) is:

Thus we have ten equations(eqs. 27.1 to 27.10) with six unknowns(a0---a5). Boundary conditions will also supply more equations. For instance two boundary conditions will add additional two equations, making a total of twelve. Following the usual procedure, we need to modify eq. 27 set so as to match the number of unknowns. Instead of merging /combining, we need to adopt a less cumbersome and easily programmable technique of defining new variables so as to match the number of equations.Thus by dividing eqs. 27.1 to 27.9 by eq. 27.10 followed by dividing eqs. 27.6 to 27.10 only with a5a5 and defining c₅ = a₅/a₅ ,c₄ = a₄/a₅, c₃ = a₃/a₅, c₂ = a₂/a₅ , c₁ = a₁/a₅,c₀ = a₀/a₅,eqs. 27 set now becomes: