Fixed Point Analytical Method for Nonlinear Differential Equations

Nonlinear problems are of great interest to physicists, mathematicians, and engineers, because most physical systems are inherently nonlinear in nature. Nonlinear problems give rise to some important phenomena, such as solitons, shock waves, chaos, and turbulence in fluid flow. Nonlinear problems are more difficult to solve than linear ones. Thus far, very few nonlinear problems are known to have a simple closed form solution. In most cases, we have to resort to an infinite series to express the solution of nonlinear problems.

where $UN$ is named as the $N$th-order approximation to $u$⁠.

Although the aforementioned idea is rather clear, there are some important questions to be answered:

1. (A)

   How do we choose the small parameter $ɛ$⁠, especially when the small parameter $ɛ$ does not obviously exist in a problem?

2. (B)

   Does the series (2) converge to the exact solution $u$?

3. (C)

   What is the influence of parameter $ɛ$ on the accuracy of $UN$⁠, especially when $ɛ$ becomes large?

It is hard to answer the preceding questions in detail. Generally speaking, the parameter $ɛ$ is usually a measure of the strength of the nonlinearity. The larger $ɛ$ is, the stronger the nonlinearity. Almost all perturbation techniques are based on the so-called small parameter assumption and the approximate solutions can be expressed in a series of small parameters, so that these perturbation techniques are applicable to the weakly nonlinear problems, and the accuracy of the perturbation series (4) becomes bad when $ɛ$ grows. What's more if the convergence of series (4) is not guaranteed, the addition of the higher-order correction term just makes the solution worse, so that the accuracy of the perturbation series cannot be improved, unlimited only by introducing more higher-order correction terms. Moreover, in some cases, the exact solution of the original problem is regular, however, the perturbation series solution (4) has singularity; even for the larger order of $UN$⁠, the stronger the singularity grows. To remove the singularity, more complicated techniques should be used. Furthermore, the perturbation series (4) usually have limited regions of validity and break down in other regions, called regions of nouniformity. To render these expansions uniformly valid, many complex and subtle techniques [1–3] should be used.

The series $ys-e$ is divergent, because the zeroth-order solution of the straightforward expansion is singular at $x=0$⁠, and the singularity grows stronger for the higher order.

which contains the so-called mixed secular term $tsint$ that tends to infinity as $t→∞$⁠. Eq. (11) is valid only for times such that $ɛt<O(1)$⁠. In the higher order approximation, the mixed secular term always exists. Moreover, $Us-e,1$ is nonperiodic, although the Duffing equation (8) denotes a cyclic motion. The breakdown in the straightforward expansion is due to its failure to account for the dependence of the frequency on the nonlinearity.

The aforementioned two examples show that the straightforward expansion is invalid for many cases and the complicated and subtle perturbation techniques, such as the PLK method and the multiple scales method [1–3], should be introduced to remove the preceding singularity, mixed secular term, and so on.

It is clear that the so-called small parameter assumption brings out all of these limitations. Therefore, it is necessary to develop a new kind of analytical method which does not depend on any small parameters at all.

where $T[x]$ is a map. The iteration procedure starts with some arbitrary initial value $x0$⁠, and the solution sequence ${xn|n=0,1,2,3,…}$ will usually converge to $x*$⁠, provided that this initial guess $x0$ is close enough to the unknown zero point $x*$⁠. In functional analysis, the zero point is named as the fixed point of the contractive map $T[x]=x-f(x)/f'(x)$⁠.

The optimal value of the relaxation factor $β$ usually is dependent on the problem to be solved.

$u*$ is called as a fixed point of the contractive map $T[u]=u-β·LC-1[A[u]]$⁠.

The fixed analytical method, which is different from the perturbation technique, is independent of any small parameters at all. Moreover, this method provides us with a convenient way to obtain a solution sequence ${un|n=0,1,2,3,…}$⁠, which can approach the exact solution of the original nonlinear equation as accurately as possible. Furthermore, the convergence of this solution sequence is so rapid that usually the first few items of the solution sequence ${un|n=0,1,2,3,…}$ can give an accurate enough approximation, which will be shown in Sec. 3.

The common solution of the linear equation $LC[u]=0$ is named as the kernel of the linear characteristic operator $LC[·]$⁠. The determination of the linear characteristic operator $LC[·]$ is not done in a systematical way, however, $LC[·]$ should obey some fundamental rules:

1. (A)

   The operator $LC[·]$ is a linear continuous bijective operator.

2. (B)

   The operator $LC[·]$ should possess as many similar properties of $A[·]$ as possible.

3. (C)

   The operator $LC[·]$ should ensure that the map $T[u]=u-β·LC-1[A[u]]$ is a contractive map.

Usually, the contraction property of the map $T[u]=u- β·LC-1[A[u]]$ is difficult to prove directly, however, it can be heuristically ensured by plotting the so called $β$-curves. It will be discussed in detail in Sec. 3.2.

In general, the steps of the fixed point analytical method applied to solve the differential equation are:

1. (A)

   To analyze the property of the original differential equation, such as the distribution of the singularity points, the asymptotic property, the continuity property, the distribution of the maximum and the minimum.

2. (B)

   To choose a basis functions system, i.e., ${ek|k=0,1,2,3,…}$⁠, which possess the preceding properties as much as possible, and construct a linear characteristic operator $LC[·]$⁠, whose kernels are members of the basis function system ${ek|k=0,1,2,3,…}$⁠. Usually, the basis function system is not unique, so we have a great amount of freedom to construct the linear characteristic operator $LC[·]$⁠.

3. (C)

   To solve the iteration procedure (15) and obtain the solution sequence ${un|n=0,1,2,3,…}$⁠, whose members should be expressed by a sum of the basis functions, i.e., $un=∑kakek,n=0,1,2,3,…$⁠, where $ak$ is the coefficient.

4. (D)

   To plot the so-called $β$-curves to decide the optimal value $β$⁠, which can greatly improve the convergence and stability of the iteration procedure.

In the same manner as Newton’s method, usually a good initial guess $u0$ will accelerate the convergence of the iteration procedure, but the question is, how to choose a good initial guess? In the framework of the FPM, any order approximation $un$ should be expressed by a sum of the basis functions, so that the initial guess of the solution is largely decided by the linear characteristic operator, the basis function system, and the initial/boundary condition. A simple and convenient initial guess $u0$ can be decided, as in the equation (17).

In this section, the fixed point analytical method is used to investigate the previous examples. All of the calculations are implemented on a laptop PC with 2 GB RAM and Intel(R) Core(TM)2 Duo 1.80GHz CPU.

where $λ$ is a free scale parameter and is decided later. The kernel of the linear characteristic equation $yker=1/(x+λ)$ satisfies the equation $LC[C/(x+λ)]=0$⁠, where $C$ is a constant of integration. It is clear that the kernel $yker$ belongs to the basis function system ${(x+λ)-k|k=0,1,2,…}$⁠. To satisfy $LC[y0]=0$ and $y0(1)=2$⁠, the initial guess $y0$ takes the form $y0=2(1+λ)/(x+λ)$⁠.

where the coefficient $an,k$ is dependent on the value of $ɛ$ and $λ$⁠.

Let us consider some typical values of the parameter $ɛ$⁠.

First, we regard the value of $λ$ as a free parameter and consider the convergence of $yn(0)$⁠, which is dependent on $λ$⁠. We plot the so-called $λ$-curves of $yn(0)$⁠. According to these $λ$-curves, we can straightforwardly and heuristically decide the valid region of $λ$⁠, which corresponds to the line segments almost parallel to the horizontal axis. In the neighborhood of $λ=0.045$⁠, the approximations $y2(0)$⁠, $y3(0)$⁠, and $y4(0)$ almost converge to the same value and the valid convergence regions of $λ$ are enlarged when the order $n$ increases, as shown in Fig. 1. In the following calculation, we set $λ=0.045$ when $ɛ=0.001$ and obtain the value of $yn(0)$ and $y'n(0)$⁠, as shown in Table 1. The comparison between the exact solution $ye(x)$ and the $n$th-order approximation $yn(x)$ is shown in Fig. 2. It is clear that a boundary layer exists at $x=0$ and the slope at $x=0$ is very large, such that $y'e(0)=-977.662$⁠. As shown in Fig. 3, the maximum relative error $([yn(x)-ye(x)]/ye(x))$ in the whole region $0≤x≤1$ is about 8.5%, 2.0%, and 0.5% for $y1(x)$⁠, $y2(x)$⁠, and $y3(x)$⁠, respectively, which shows that the convergence of the solution sequence ${yn|n=0,1,2,3,…}$ to $ye(x)$ is rapid and uniform.

As shown in Fig. 4, in the neighborhood of $λ=2.5$⁠, the approximations $y2(0)$⁠, $y3(0)$⁠, and $y4(0)$ almost converge to the same value and the valid convergence regions of $λ$ are enlarged when the order $n$ increases. We set $λ=2.5$ when $ɛ=1$ and obtain the value of $yn(0)$ and $y'n(0)$⁠, as shown in the Table 2.

The comparison between the exact solution $ye(x)$ and the $n$th-order approximation $yn(x)$ is shown in Fig. 5. The slope at $x=0$ is relatively small $y'e(0)=-0.591752$⁠, and the boundary layer at $x=0$ disappears. As shown in Fig. 6, the maximum relative error (⁠$[yn(x)-ye(x)]/ye(x)$⁠) in the whole region $0≤x≤1$ is about 0.7%, 0.05%, and 0.006% for $y1(x)$⁠, $y2(x)$⁠, and $y3(x)$⁠, respectively, which, again, shows that the convergence of the solution sequence ${yn|n=0,1,2,3,…}$ to $ye(x)$ is rapid and uniform.

Using a similar method, we set the scale parameter $λ=2000$ when $ɛ=1000$ and obtain the values of $yn(0)$ and $y'n(0)$⁠, as shown in Table 3. The comparison between the exact solution $ye(x)$ and the $n$th-order approximation $yn(x)$ is shown in Fig. 7. The slope at $x=0$ is even smaller $y'e(0)=-5.00125×10-4$⁠, and the profile of $y(x)$ is almost a horizontal line. As shown in Fig. 8, the maximum relative error in the whole region $0≤x≤1$ is so tiny, about $1.2×10-5%$ for $y1(x)$⁠, that $y1(x)$ can give an accurate enough approximation to $y(x)$ when $ɛ=1000$⁠.

agrees well with the exact value $ye(0)$ in the whole region $0<ɛ<∞$⁠, as shown in Fig. 9. The maximum relative error $(y2(0)/ye(0)-1)$ occurs when $ɛ=1/6$

It is shown that the free scale parameter $λ$ provides us with a convenient way to ensure the convergence of the solution sequence. Different from the perturbation approximations, Eqs. (26) and (28) are valid for all possible parameters $0<ɛ<∞$⁠, and thus, they are independent of the small parameters.

We have known that the straightforward expansion (11) has the mixed secular term, which is valid only for $ɛt<O(1)$⁠. Without loss of generality, we first consider the case $ɛ=1$⁠. The procedure for the large value of $ɛ$ is similar.

where $C1$ and $C2$ are constant. Because $uker,2=sin(τ)$ does not belong to the basis function system ${cos(kτ)|k=0,1,2,3,…}$⁠, we have to discard $uker,2=sin(τ)$ and set the constant $C2=0$⁠.

where $β$ is the relaxation factor, which is a nonzero free parameter used to improve the convergence, and $Ωn$ is the $n$th-order approximation to the exact value $Ω$⁠.