A power series method for solving ordinary and partial differentials equations motivated by domain growth

TL;DR

This paper introduces a power series approach for solving various ordinary and partial differential equations, demonstrated through multiple examples including linear and nonlinear cases, with boundary conditions and general solutions provided.

Contribution

The paper presents a novel power series method applicable to both ordinary and partial differential equations, including nonlinear and boundary value problems.

Findings

Successfully solves diverse differential equations using the power series method

Demonstrates applicability to linear and nonlinear PDEs and ODEs

Provides solutions with boundary conditions and general forms

Abstract

In this work we present a power series method for solving ordinary and partial differential equations. To demonstrate our method we solve a system of ordinary differential equations describing the movement of a random walker on a one-dimensional lattice, two nonlinear ordinary differential equations, a wave and diffusion equation (linear partial differential equations), and a nonlinear partial differential equation (quasilinear). The inclusion of boundary conditions and the general solutions to other equations of interest are included in the Supplementary material.

Equations

$$\frac{\mathrm{d}{p}_i}{\mathrm{d}t}=\left(\frac{P_m}{2}\right)\left(p_{i-1}-2p_i+p_{i+1}\right).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\sum _{n=0}^{\infty }{p}_i^n=\sum _{n=0}^{\infty }\left(\frac{P_m}{2}\right)\left(p_{i-1}^n-2p_i^n+p_{i+1}^n\right).$$

$$p_i=\sum _{n=0}^{\infty }{p}_i^n.$$

$$\frac{\mathrm{d}{p}_i^0}{\mathrm{d}t}=-\beta p_i^0,$$

$$\frac{\mathrm{d}{p}_i^n}{\mathrm{d}t}=-\beta p_i^n+\beta p_i^{n-1}+\left(\frac{P_m}{2}\right)\left(p_{i-1}^{n-1}-2p_i^{n-1}+p_{i+1}^{n-1}\right),\forall n>0.$$

$$\displaystyle\frac{\mathrm{d}p^{n}_{i}}{\mathrm{d}t}=-P_{g}p^{n}_{i}+P_{g}p^{n-1}_{i}+\left(\frac{P_{m}}{2}\right)\Big{(}p^{n}_{i-1}-2p^{n}_{i}+p^{n}_{i+1}\Big{)},$$

$$\frac{\mathrm{d}{p}_i^0}{\mathrm{d}t}=0,$$

$$\frac{\mathrm{d}{p}_i^n}{\mathrm{d}t}=\left(\frac{P_m}{2}\right)\left(p_{i-1}^{n-1}-2p_i^{n-1}+p_{i+1}^{n-1}\right),\forall n>0.$$

$$\frac{\mathrm{d}{p}_i^n}{\mathrm{d}t}=-\beta p_i^n+(\beta -P_m)p_i^{n-1}+\left(\frac{P_m}{2}\right)\left(p_{i-1}^{n-1}+p_{i+1}^{n-1}\right),\forall n>0.$$

$$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\sum^{\infty}_{n=0}p^{n}_{i}=\sum^{\infty}_{n=0}\left(\frac{P_{m}}{2}\right)\Big{(}p^{n}_{i-1}-2p^{n}_{i}+p^{n}_{i+1}\Big{)}.$$

$$\frac{\mathrm{d}{p}_i}{\mathrm{d}t}=\left(\frac{P_m}{2}\right)\left(p_{i-1}-2p_i+p_{i+1}\right),$$

$$\frac{\mathrm{d}{p}_i^0}{\mathrm{d}t}=-\beta p_i^0.$$

$$p_i^0=A_i\exp (-\beta t),$$

$$\frac{\mathrm{d}{p}_i^1}{\mathrm{d}t}=-\beta p_i^1+(\beta -P_m)A_i\exp (-\beta t)+\left(\frac{P_m}{2}\right)\left(A_{i-1}+A_{i+1}\right)\exp (-\beta t).$$

$$p_i^n(t)=\exp (-\beta t)\left(\frac{t^n}{n!}\right)\left[\sum _{j=0}^n(\beta -P_m{)}^{n-j}(P_m{)}^j\frac{1}{2^j}\left(\genfrac{}{}{0ex}{}{n}{j}\right)\left(\sum _{k=-j}^{-j:2:j}\left(\genfrac{}{}{0ex}{}{j}{\frac{(k+j)}{2}}\right)A_{i+k}\right)\right],\forall n\ge 0.$$

$$\sum _{n=0}^{\infty }{p}_i^n(t)=\sum _{n=0}^{\infty }\exp (-\beta t)\left(\frac{t^n}{n!}\right)\left[\sum _{j=0}^n(\beta -P_m{)}^{n-j}(P_m{)}^j\frac{1}{2^j}\left(\genfrac{}{}{0ex}{}{n}{j}\right)\left(\sum _{k=-j}^{-j:2:j}\left(\genfrac{}{}{0ex}{}{j}{\frac{(k+j)}{2}}\right)A_{i+k}\right)\right],$$

$$\frac{\mathrm{d}{p}_1^n}{\mathrm{d}t}$$

$$\text{mathmakebox}[\text{widthof}=][c]⋮$$

$$\frac{\mathrm{d}{p}_i^n}{\mathrm{d}t}$$

$$\text{mathmakebox}[\text{widthof}=][c]⋮$$

$$\frac{\mathrm{d}{p}_N^n}{\mathrm{d}t}$$

$$p_i^n(t)=\left(\frac{t}{n}\right)\left[(\beta -P_m)p_i^{n-1}(t)+\left(\frac{P_m}{2}\right)(p_{i+1}^{n-1}(t)+p_{i-1}^{n-1}(t))\right],\forall n>0,$$

$$p_i^0(t)=A_i\exp (-\beta t).$$

$$\sum _{n=0}^{\infty }{p}_i^n(t)=p_i^0(t)+\sum _{n=1}^{\infty }\left(\frac{t}{n}\right)\left[(\beta -P_m)p_i^{n-1}(t)+\left(\frac{P_m}{2}\right)(p_{i+1}^{n-1}(t)+p_{i-1}^{n-1}(t))\right].$$

$$p_1^n(t)=\left(\frac{t}{n}\right)\left[\left(\beta -\frac{P_m}{2}\right)p_1^{n-1}(t)+\left(\frac{P_m}{2}\right)(p_2^{n-1}(t))\right]\forall n>0,$$

$$p_N^n(t)=\left(\frac{t}{n}\right)\left[\left(\beta -\frac{P_m}{2}\right)p_N^{n-1}(t)+\left(\frac{P_m}{2}\right)(p_{N-1}^{n-1}(t))\right]\forall n>0.$$

$$\sum _{n=0}^{\infty }{p}_1^n(t)=p_1^0(t)+\sum _{n=1}^{\infty }\left(\frac{t}{n}\right)\left[\left(\beta -\frac{P_m}{2}\right)p_1^{n-1}(t)+\left(\frac{P_m}{2}\right)(p_2^{n-1}(t))\right],$$

$$\sum _{n=0}^{\infty }{p}_N^n(t)=p_N^0(t)+\sum _{n=1}^{\infty }\left(\frac{t}{n}\right)\left[\left(\beta -\frac{P_m}{2}\right)p_N^{n-1}(t)+\left(\frac{P_m}{2}\right)(p_{N-1}^{n-1}(t))\right],$$

$$\frac{\mathrm{d}p}{\mathrm{d}t}=\gamma p(1-p),$$

$$p(t)=\frac{C_1\exp (\gamma t)}{1-C_1+C_1\exp (\gamma t)},$$

$$\frac{\mathrm{d}{p}^n(t)}{\mathrm{d}t}=-\beta p^n(t)+\beta p^{n-1}(t)+\gamma p^{n-1}(t)-\gamma p^{n-1}(t)\sum _{n=0}^{\infty }{p}^n(t),\forall n>0,$$

$$\frac{\mathrm{d}{p}^0(t)}{\mathrm{d}t}=-\beta p^0(t).$$

Full text

Keywords: Differential equations, ordinary, partial, power series solutions, growth.

1 Introduction

We present a method that can be used to generate power series solutions for linear and nonlinear ordinary differential equations (ODEs) and partial differential equations (PDEs). The method we present relies on a separation of variables in a system of equations we construct, and generates a power series weighted by coefficients written in terms of the initial condition. This method was motivated by work conducted on the mathematical modeling of domain growth [1].

The outline of this work is as follows: In Section 2.1 we demonstrate how our method can be used to generate power series solutions for linear and nonlinear ODEs. In Sections 2.2 and 2.3 we demonstrate how this method can be extended to linear and nonlinear PDEs, including the implementation of boundary conditions in the linear PDE case. We finish with a short discussion of the method presented in this work in Section 3.

2 Results

2.1 Solving a system of ordinary differential equations describing the movement of a random walker on a one-dimensional lattice

A one-dimensional lattice with periodic boundaries (a ring) is displayed in Fig. 1.

The following equation describes the time evolution of the probability that an unbiased excluding random walker occupies site $i$ on a one-dimensional periodic lattice111A derivation of Eq. (1) can be found in the Supplementary material (SM1).:

[TABLE]

In Eq. (1) $p_{i}$ is the probability a random walker is situated at site $i$ at time $t$, and $P_{m}$ is the rate at which the random walker attempts to move to an adjacent site on the lattice.

We now write Eq. (1) in the following manner

[TABLE]

That is, we postulate that $p_{i}$ can be written as the infinite series

[TABLE]

We now decompose Eq. (2) into the following infinite system of equations

[TABLE]

and

[TABLE]

Notice in Eq. (5) that we have written the terms associated with the movement of the random walker in terms of $n-1$, not $n$.222It has previously been shown [1] that if the lattice is growing, Eq. (5) would be written as

\displaystyle\frac{\mathrm{d}p^{n}_{i}}{\mathrm{d}t}=-P_{g}p^{n}_{i}+P_{g}p^{n-1}_{i}+\left(\frac{P_{m}}{2}\right)\Big{(}p^{n}_{i-1}-2p^{n}_{i}+p^{n}_{i+1}\Big{)},

where the motility terms have the same $n$ as the time derivative, and $P_{g}$ is the rate at which the lattice grows. This observation is what initially motivated the work we present here. We include the parameter $\beta$ in Eqs. (4) and (5) as a ‘shape’ parameter, however, its inclusion in Eqs. (4) and (5) is not necessary and $\beta$ can be set to zero if desired.333If $\beta=0$ then Eqs. (4) and (5) would be

$\displaystyle\frac{\mathrm{d}p^{0}_{i}}{\mathrm{d}t}=0,$

(6)

and

$\displaystyle\frac{\mathrm{d}p^{n}_{i}}{\mathrm{d}t}=\left(\frac{P_{m}}{2}\right)\left(p^{n-1}_{i-1}-2p^{n-1}_{i}+p^{n-1}_{i+1}\right),\ \ \ \forall\ n>0.$

(7) Finally, we simplify Eq. (5) to obtain

[TABLE]

It is readily apparent that if we sum Eq. (4) and Eq. (5) for all $n>0$ we obtain

[TABLE]

This means if we substitute Eq. (3) into Eq. (9) we arrive at

[TABLE]

which recapitulates Eq. (1).

The decomposition of Eq. (1) into the infinite system of equations contained in Eqs. (4) and (5) is straightforward to solve. To see this consider the initial equation, Eq. (4),

[TABLE]

As the initial equation does not ‘inherit’ any terms associated with the movement of the random walker, Eq. (11) admits the simple solution:

[TABLE]

where $A_{i}$ is the initial value of site $i$.

Equation (12) can be placed in Eq. (8) for $n=1$, so that Eq. (8) becomes

[TABLE]

This means Eq. (13) is now also straightforward to solve. Repeated application of this process admits the following recurrence formula as a solution for Eq. (8):

[TABLE]

Therefore, the probability of site $i$ being occupied by a walker at time $t$ is given by

[TABLE]

in accordance with Eq. (3).

In Fig. 2 we demonstrate that the solution given by Eq. (15) matches the evolution of the ensemble average of the discrete model excellently at our chosen level of truncation (for the algorithm used in the discrete model see the Supplementary material (SM2)). In Fig. 3 we display what we refer to as ‘streams’, given by Eq. (14), for a single value of the shape parameter $\beta$. In the Supplementary material (SM3) we display streams for different values of $\beta$, which demonstrates how $\beta$ influences the shape of the streams that compose the solution given by Eq. (15).

2.1.1 Boundary conditions

It is possible to implement boundary conditions in equations describing the movement of a random walker on a one-dimensional lattice with the method we are presenting. For an unbiased random walker on a one-dimensional lattice with no-flux boundary conditions the equations describing the probability of finding a walker at a given site are:

[TABLE]

Following the same procedure as we did for Eqs. (14) and (15) we find the following recurrence relation for sites not situated on the boundary

[TABLE]

with

[TABLE]

The reader will notice that we have written Eq. (17) in a more economical form than Eq. (14). It is possible to write Eq. (17) in the same manner as Eq. (14), and if done so each site on the lattice will have its own recurrence formula describing the probability of a random walker being located at that site at time $t$.

From Eq. (17) the evolution with respect to time of the probability of site $i$ being occupied by a walker is given by

[TABLE]

The recurrence relations for sites situated on the boundary are

[TABLE]

and

[TABLE]

Therefore, the evolution with respect to time of the probability of sites $1$ and $N$ being occupied by a walker are given by

[TABLE]

and

[TABLE]

respectively.

In Fig. 4 the solutions given by Eqs. (19), (22) and (23) are displayed. It can be seen that Eqs. (19), (22) and (23) and the ensemble average from the discrete model match excellently. The algorithm for the discrete model can be found in the Supplementary material (SM2).

We provide a final example of our method being applied to a linear ODE in the Supplementary material (SM4).

2.1.2 Nonlinear ordinary differential equations

We now apply our method to nonlinear ODEs. This allows us to demonstrate another important aspect of our method. Initially we solve

[TABLE]

where $\gamma>0$. The analytic solution of Eq. (24) is

[TABLE]

where $C_{1}$ is the value of $p(t)$ at $t=0$.

To solve Eq. (24) in our framework we proceed as follows. To begin with we decompose Eq. (24) into

[TABLE]

with

[TABLE]

It is evident that Eq. (26) cannot be solved in the same iterative manner as Eq. (5) due to the due to the nonlinear term present on its right-hand-side.444One might think the Eq. (26) should be written as

$\displaystyle\frac{\mathrm{d}p^{n}(t)}{\mathrm{d}t}=-\beta p^{n}(t)+\beta p^{n-1}(t)+\gamma p^{n-1}(t)-\gamma p^{n-1}(t)p^{n-1}(t),\ \ \ \forall\ n>0,$

(28)

but this is incorrect as each stream needs to be multiplied by every other stream to account for the nonlinearity in Eq. (24). To circumvent this we sum Eq. (26) for all $n\geq 0$ to obtain

[TABLE]

and then decompose Eq. (29) in the following manner:

[TABLE]

[TABLE]

[TABLE]

so that in general

[TABLE]

The decomposition of Eq. (29) into the equations contained in Eqs. (30) and (33) allows us to solve the unknowns iteratively, and it is straightforward to demonstrate that summing Eqs. (30) and (33) for all $n>0$ returns Eq. (29). In Fig. 5 (a) we compare the solution of Eq. (33) with the analytical solution Eq. (25).

It is also possible to derive power series solutions in terms of simple functions for nonlinear ODEs with the method we are presenting. For instance, the nonlinear ODE

[TABLE]

has the following power series solution

[TABLE]

where

[TABLE]

In Fig. 5 (b) we compare Eq. (35) with the analytical solution of Eq. (34),

[TABLE]

The details of how to derive Eq. (35) are given in the Supplementary material (SM5).

2.2 Solving a linear partial differential equation

We now demonstrate that the method we are presenting is also applicable to PDEs. We start by applying this method to linear PDEs, for instance the wave equation:

[TABLE]

Motivated by the previous section we write Eq. (38) as

[TABLE]

with

[TABLE]

As before we initially solve Eq. (40),

[TABLE]

It can be seen in Eq. (41) that our method relies on the separation of spatial and temporal variable in the initial equation. For general $u^{n}$ we obtain

[TABLE]

If $\beta=0$ we have

[TABLE]

and each stream is a polynomial in $t$ weighted by coefficients written in terms of the initial condition555In the case of Eq. (43) it is evident we have simply derived a Taylor series expansion [2].. From Eq. (42) the general solution to Eq. (38) is

[TABLE]

A simple initial condition for Eq. (44) is

[TABLE]

and if we substitute Eq. (45) into Eq. (42) we obtain

[TABLE]

It is also straightforward to solve the diffusion equation, which is

[TABLE]

Solving Eq. (47) in a similar manner to how we solved Eq. (38) we obtain

[TABLE]

If $\beta=0$ Eq. (48) is

[TABLE]

Therefore, from Eq. (48) our general solution to Eq. (47) is

[TABLE]

If we use $A(x)$ = sin($x$) for the initial condition in Eq. (48) we obtain

[TABLE]

This method is trivially extendable to two-dimensional linear PDEs, the details of which are given in the Supplementary material (SM6). It is also possible to implement boundary conditions, and this is also demonstrated in the Supplementary material (SM7).

2.3 Solving a nonlinear partial differential equation

Finally, we demonstrate that this method is also extendable to nonlinear PDEs. For instance, the quasilinear inviscid Burgers equation. The inviscid Burgers equation is

[TABLE]

where $\alpha$ is a constant. We begin by writing Eq. (52) in the following manner

[TABLE]

with

[TABLE]

As in the case of nonlinear ODEs we have to multiply the $n^{th}$ stream by all other streams (including itself) to account for the nonlinearity in Eq. (52). We then decompose Eqs. (53) and (54) in the following manner:

[TABLE]

with

[TABLE]

and

[TABLE]

In Fig. 6 the solution of Eqs. (55)-(57) is compared with the solution of Eq. (52) before the onset of the multivalue behaviour that the solution of Eq. (52) exhibits. We use symbolic integration in Matlab to compute Eqs. (55)-(57). It should be readily apparent how to extend this method to more complicated nonlinear PDEs.

3 Discussion

We have presented a power series method for solving both linear and nonlinear ODEs and PDEs. We finish by detailing some issues with the method we have introduced in this work.

Our main criticism of the work we have presented is that in the case of some nonlinear equations presented in this work we have not supplied solutions for the $n^{th}$ stream written in terms of simple functions. For instance, Eqs. (24) and (52). The method presented here would be most useful if an efficient means of writing the power series solutions for nonlinear equations became evident, which would allow analysis to be directly carried out on these solutions. The Supplementary material (SM5) shows that in some cases of nonlinear equations it is possible to write the $n^{th}$ stream of its solution in terms of simple functions, however, a way to generalise the approach used on Eq. (34) has not yet become apparent to the authors.

It is also important to acknowledge that we have not dealt with the issue of convergence in the power series we have presented. It is obvious to say that the convergence of these power series, and their radius of convergence, will depend on the initial conditions of the equation, and the equation itself [2]. However, a more general treatment on the convergence of the methods presented here is certainly required. Finally, a word on the role of the shape parameter $\beta$. Its role may seem somewhat superfluos, however, it is a simple way to circumvent numerical issues when the value of streams that compose solutions becomes too large for a standard computer to accurately represent. It also means that the value of the streams composing a solution can be made positive for a given interval of interest by selecting the appropriate value of $\beta$, and so provides another analytic tool to utilise when employing the methods presented here.

Acknowledgements

RJHR would like to thank Kit Yates, Ruth Baker and Pierre Boutillier for helpful discussions. The author declares no competing interests.

Supplementary material

SM1: The derivation of Equation (1) in the main text.

We derive Eq. (1) in the following manner. The probability that an unbiased excluding random walker occupies site $i$ on a one-dimensional periodic lattice at time $t+\delta t$ is given by

[TABLE]

In Eq. (58) $p_{i-1,i}(A,0;t)$ is the second-order probability that site $i-1$ and $i$ are occupied and unoccupied, respectively, at time $t$. The other second-order terms in Eq. (58) have similar meanings. If we rearrange Eq. (58) and take $\delta t\rightarrow 0$ in the limit we obtain

[TABLE]

We now remove the second-order terms in Eq. (59) by making the following closure

[TABLE]

If we place Eq. (60) in Eq. (61) we obtain

[TABLE]

If we drop the explicit ‘$A$’ and ‘$t$’ from our notation in Eq. (61) we recapitulate Eq. (1).

SM2: Algorithm for discrete random-walk

We use a discrete random-walk model on a one-dimensional regular lattice with lattice spacing $\Delta$ [3] and length $N$, where $N$ is an integer describing the number of lattice sites. Simulations are performed with either periodic boundary or no-flux conditions. Each random walker is assigned to a lattice site, from which it can move into an adjacent site. If an agent attempts to move into a site that is already occupied, the movement event is aborted. This process, whereby only one agent is allowed per site, is generally known as an exclusion process. Time is evolved continuously, and random walker movements are attempted in accordance with the Gillespie algorithm [4]. Attempted agent movement events occur with rate $P_{m}$ per unit time. The initial conditions of the discrete model are provided in the main text when necessary.

SM3: The effect of different values of $\beta$ in Eq. (15) in the main text.

In Fig. 7 we display streams, Eq. (14), for different values of $\beta$. This demonstrates how $\beta$ influences the shape of the streams that compose the solution given by Eq. (15).

SM4: Linear ordinary differential equation

We provide the solution to the following linear ODE

[TABLE]

which is linear in $q$. The analytic solution of Eq. (62) is

[TABLE]

where $C_{2}$ is the value of $q(t)$ at $t=0$. To solve Eq. (62) in our framework we rewrite it as

[TABLE]

with

[TABLE]

In Fig. 6 (a) the solution of Eqs. (64) and (65) is compared with the analytical solution Eq. (63).

SM5: Nonlinear ordinary differential equation

To solve Eq. (34) we proceed in the following manner. We begin with

[TABLE]

and

[TABLE]

and

[TABLE]

In this derivation we assume $\beta=0$ for simplicity. Initially, we rewrite Eq. (66) as

[TABLE]

This means

[TABLE]

which gives

[TABLE]

Therefore

[TABLE]

because $c_{1}=-y^{0}\operatorname{log}(y^{0})$. If we place Eq. (72) in Eq. (68) we obtain

[TABLE]

which we can integrate to obtain

[TABLE]

Now we recognise

[TABLE]

If we integrate Eq. (75), and then solve for $y^{3}$ in a similar manner we obtain the following power series solution for $y$

[TABLE]

Alternatively, we can place Eq. (74) into Eq. (72) to obtain

[TABLE]

which recapitulates the analytical solution to Eq. (66), given as Eq. (37).

SM6: Two-dimensional linear partial differential equation

The general solution for the two-dimensional linear diffusion equation $(\beta=0)$ is

[TABLE]

SM7: Boundary conditions for linear partial differential equation

We now demonstrate how to implement boundary conditions in linear PDEs with our method. A simple example is for the diffusion equation

[TABLE]

with

[TABLE]

and

[TABLE]

If we implemented Neumann boundary conditions these would take the form

[TABLE]

To implement our boundary conditions we proceed as follows: As we already have the general solution for Eq. (48) we can take its partial derivative with respect to time to obtain666A simpler way to obtain the solution for the given initial condition is the following. Implementing initial condition in Eq. (79) gives:

$\displaystyle\frac{\partial u^{L+\delta L}}{\partial t}=-\beta u^{L+\delta L}+\beta u^{L},$

(83)

with the (straightforward) solution being

$\displaystyle\sum^{\infty}_{i=0}u^{L+i\delta L}(x,t)=\sum^{\infty}_{i=0}A(x)\frac{(\beta t)^{i}}{i!}\operatorname{exp}(-\beta t)=A(x)\sum^{\infty}_{i=0}\frac{(\beta t)^{i}}{i!}\operatorname{exp}(-\beta t)=A(x),$

(84)

as the Poisson distribution sums to identity.

[TABLE]

which means

[TABLE]

Integrating Eq. (86) with respect to $x$ gives

[TABLE]

and integrating Eq. (87) with respect to $x$ gives

[TABLE]

If we apply the boundary conditions and initial condition to Eq. (88) we obtain

[TABLE]

as the solution to Eq. (79). In this solution we define

[TABLE]

when $n=0$ to avoid a singularity when $t=0$.