On the generalized method of lines and its proximal explicit and hyper-finite difference approaches

TL;DR

This paper introduces proximal explicit and hyper-finite difference methods for the generalized method of lines, aiming to reduce solution error in PDE discretization through novel domain decomposition and boundary condition techniques.

Contribution

It develops new proximal and hyper-finite difference approaches for the generalized method of lines, enhancing accuracy and computational efficiency in PDE solutions.

Findings

Proximal explicit approach reduces solution error for small parameters.

Hyper-finite differences enable domain decomposition without increasing error.

Numerical examples demonstrate method effectiveness.

Abstract

This article firstly develops a proximal explicit approach for the generalized method of lines. In such a method, the domain of the PDE in question is discretized in lines and the equation solution is written on these lines as functions of the boundary conditions and domain shape. The main objective of introducing a proximal formulation is to minimize the solution error as a typical parameter $\varepsilon>0$ is too small. In a second step we present another procedure to minimize this same error, namely, the hyper-finite differences approach. In this last method the domain is divided in sub-domains on which the solution is obtained through the generalized method of lines allowing the parameter $\varepsilon>0$ to be very small without increasing the solution error. The solutions for the sub-domains are connected through the boundary conditions and the solution of the partial differential equation in question on the node lines which separate the sub-domains. In the last sections of each text part we present the concerning softwares and perform numerical examples.

Equations

$$\epsilon {\nabla }^{2}u+g(u)+f=0,\text{ in }\Omega \subset {\mathbb{R}}^2,$$

$$u=0\text{ on }{\Gamma }_0\text{ and }u=u_f,\text{ on }{\Gamma }_1.$$

$${\Gamma }_1=2{\Gamma }_0,$$

$${\Gamma }_0=\{(\theta ,r(\theta ))|0\le \theta \le 2\pi \},$$

$$\frac{{\partial }^{2}u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}+\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\theta }^2}+g(u)+f=0,\text{ in }\Omega ,$$

$$u=0\text{ on }{\Gamma }_0\text{ and }u=u_f,\text{ on }{\Gamma }_1.$$

$$t=\frac{r}{r(\theta )}.$$

$$u(r,\theta )=\bar{u}(t,\theta ),$$

$$\frac{{\partial }^{2}u}{\partial t^2}$$

$$f_1(\theta )=\frac{-r^{\prime }(\theta )}{r(\theta )},$$

$$f_2(\theta )=1+\frac{f_1^{\prime }(\theta )}{1+f_1(\theta {)}^2},$$

$$f_3(\theta )=\frac{2f_1(\theta )}{1+f_1(\theta {)}^2},$$

$$f_4(\theta )=\frac{1}{1+f_1(\theta {)}^2}.$$

$$\frac{u_{n+1}-2u_n+u_{n-1}}{d^2}+\frac{(u_n-u_{n-1})}{d}\frac{1}{t_n}{f}_2(\theta )$$

$$+\frac{\partial (u_n-u_{n-1})}{\partial \theta }\frac{1}{t_{n}d}{f}_3(\theta )+\frac{{\partial }^2{u}_n}{\partial {\theta }^2}\frac{f_4(\theta )}{t_n^2}$$

$$+f_5(\theta )\left(g(u_n)\frac{1}{\epsilon }+f_n\frac{1}{\epsilon }\right)=0,$$

$$u_n=T_n(u_{n-1},u_n,u_{n+1}),$$

$$T_n(u_{n-1},u_n,u_{n+1})$$

$$\Vert T(x_1)-T(x_2){\Vert }_U\le \alpha \Vert x_1-x_2{\Vert }_U,\forall x_1,x_2\in C.$$

$$\Vert T(x_1)-T(x_2){\Vert }_U\le \underset{x\in C}{\sup }\{\Vert T^{\prime }(x)\Vert \}\Vert x_1-x_2{\Vert }_U,\forall x_1,x_2\in C.$$

$$x_1=T(x_0)\text{ and }{x}_{k+1}=T(x_k),\forall k\in \mathbb{N}$$

$$x_k\to \tilde{x},\text{in norm, as }k\to +\infty .$$

$$u_{n+1}-2u_n+u_{n-1}+\frac{(u_n-u_{n-1})}{d}\frac{1}{t_n}{f}_2(\theta )d^2$$

$$T_n(u_n,u_{n+1},u_{n-1})$$

$$u_n=T_n(u_{n-1},u_n,u_{n+1}),$$

$$u_1=T_1(0,u_1,u_2).$$

$$x_{k+1}=T_1(0,x_k,u_2),$$

$$u_1=\underset{k\to \infty }{\lim }{x}_k=g_1(u_2).$$

$$u_1=g_1(u_2).$$

$$u_2=T_2(u_1,u_2,u_3),$$

$$u_2=T_2(g_1(u_2),u_2,u_3).$$

Full text

1 Introduction

This article develops two improvements relating the generalized method of lines. In our previous publications [2, 4], we highlight the method there addressed may present a relevant error as a parameter $\varepsilon>0$ is too small, that is, as $\varepsilon$ is about 0.01, 0.001 or even smaller.

In the present section we develop a solution for such a problem through a proximal formulation suitable for a large class of non-linear elliptic PDEs.

At this point we reintroduce the generalized method of lines, originally presented in F.Botelho [2]. In the present context we add new theoretical and applied results to the original presentation. Specially the computations are all completely new. Consider first the equation

[TABLE]

with the boundary conditions

[TABLE]

From now on we assume that $u_{f}$, $g$ and $f$ are smooth functions (we mean $C^{\infty}$ functions), unless otherwise specified. Here $\Gamma_{0}$ denotes the internal boundary of $\Omega$ and $\Gamma_{1}$ the external one. Consider the simpler case where

[TABLE]

and suppose there exists $r(\theta)$, a smooth function such that

[TABLE]

being $r(0)=r(2\pi)$.

In polar coordinates the above equation may be written as

[TABLE]

and

[TABLE]

Define the variable $t$ by

[TABLE]

Also defining $\bar{u}$ by

[TABLE]

dropping the bar in $\bar{u}$, equation (1) is equivalent to

[TABLE]

in $\Omega$. Here $f_{2}(\theta),\;f_{3}(\theta)$, $f_{4}(\theta)$ and $f_{5}(\theta)$ are known functions.

More specifically, denoting

[TABLE]

we have:

[TABLE]

[TABLE]

and

[TABLE]

Observe that $t\in[1,2]$ in $\Omega$. Discretizing in $t$ (N equal pieces which will generate N lines ) we obtain the equation

[TABLE]

$\forall n\in\{1,...,N-1\}$. Here, $u_{n}(\theta)$ corresponds to the solution on the line $n$. Thus we may write

[TABLE]

where

[TABLE]

1.1 Some preliminaries results and the main algorithm

Now we recall a classical definition.

Definition 1.1**.**

Let $C$ be a subset of a Banach space $U$ and let $T:C\rightarrow C$ be an operator. Thus $T$ is said to be a contraction mapping if there exists $0\leq\alpha<1$ such that

[TABLE]

Remark 1.2**.**

Observe that if $\|T^{\prime}(x)\|_{U}\leq\alpha<1,$ on a convex set $C$ then $T$ is a contraction mapping, since by the mean value inequality,

[TABLE]

The next result is the base of our generalized method of lines. For a proof see [3].

Theorem 1.3** (Contraction Mapping Theorem).**

Let $C$ be a closed subset of a Banach space $U$. Assume $T$ is contraction mapping on $C$, then there exists a unique $\tilde{x}\in C$ such that $\tilde{x}=T(\tilde{x})$. Moreover, for an arbitrary $x_{0}\in C$ defining the sequence

[TABLE]

we have

[TABLE]

To obtain a fixed point for each $T_{n}$ indicated in (5) is perfectly possible if $\varepsilon\approx\mathcal{O}(1).$ However if $\varepsilon>0$ is small, the error in this process may be relevant.

To solve this problem, firstly we propose the following algorithm,

Choose $K\approx 30-80$ and set $u_{0}=\mathbf{0}.$ 2. 2.

Calculate $u=\{u_{n}\}$ by solving the equation

[TABLE]

Such an equation is solved through the Banach fixed point theorem, that is, defining

[TABLE]

equation (14) stands for

[TABLE]

so that for $n=1$ we have

[TABLE]

We may use the Contraction Mapping Theorem to calculate $u_{1}$ as a function of $u_{2}$. The procedure would be,

1. (a)

set $x_{1}=u_{2},$ 2. (b)

obtain recursively

[TABLE] 3. (c)

and finally get

[TABLE]

Thus, we have obtained

[TABLE]

We can repeat the process for $n=2$, that is, we can solve the equation

[TABLE]

which from above stands for

[TABLE]

The procedure would be:

1. (a)

Set $x_{1}=u_{3}$, 2. (b)

calculate

[TABLE] 3. (c)

obtain

[TABLE]

We proceed in this fashion until obtaining

[TABLE]

Being $u_{f}$ known we have obtained $u_{N-1}$ . We may then calculate

[TABLE]

[TABLE]

and so on, up to finding

[TABLE]

Thus this part of the problem is solved. 3. 3.

Set $u_{0}=u$ and go to item 2 up to the satisfaction of an appropriate convergence criterion.

Remark 1.4**.**

Here we consider some points concerning the convergence of the method.

In the next lines the norm indicated refers to the infinity one for $C([0,2\pi];\mathbb{R}^{N-1})$. In particular for $n=1$ from above we have:

[TABLE]

that is

[TABLE]

Hence, denoting

[TABLE]

and

[TABLE]

for $N$ sufficiently big we may obtain

[TABLE]

and by induction

[TABLE]

so that we would have

[TABLE]

This last calculation is just to clarify that the procedure of obtaining the relation between consecutive lines through the contraction mapping theorem is well defined.

1.2 A numerical example, the proximal explicit approach

In this section we present a numerical example. Consider the equation

[TABLE]

where, for a Ginburg-Landau type equation (see [1, 5] for the corresponding models in physics),

[TABLE]

with the boundary conditions

[TABLE]

where $\Omega=\{(r,\theta)\;:\;1\leq r\leq 2,\;0\leq\theta\leq 2\pi\}$,

[TABLE]

[TABLE]

Through the generalized method of lines, for $N=10$ (10 lines), $d=1/N$ in polar coordinates and finite differences (please see [6] for general schemes in finite differences), equation (17), stands for

[TABLE]

$\forall n\in\{1,\ldots,N-1\}.$

At this point we present, through the generalized method of lines, the concerning algorithm which may be for the softwares maple or mathematica.

In this software, $x$ stands for $\theta$.

[TABLE]

[TABLE]

At this point we present the expressions for 10 lines, firstly for $\varepsilon=1$ and $K=0.$ In the next lines $x$ stands for $\theta$.

For each line $u[n]$ we have obtained,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In the next lines we present the results relating the software indicated, with $\varepsilon=0.01$ and $K=70.$

For each line $u[n]$ we have obtained,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Remark 1.5**.**

Observe that since $\varepsilon=0.01$ the solution is close to the constant value $1.3247$ along the domain, which is an approximate solution of equation $-u^{3}+u+1.0=0.$

2 The hyper-finite differences approach

In the last sections we have introduced a method to minimize the solution error for a large class of PDEs as a typical parameter $\varepsilon>0$ is small. The main idea there presented it consisted of a proximal formulation combined with the generalized method of lines.

In the present section we develop another new solution for such a problem also for a large class of non-linear elliptic PDEs, namely, the hyper-finite differences approach. We believe the result here developed is better than those of the previous sections. Indeed in the present method, in general the convergence is obtained more easily.

The idea here is to divide the domain interval, concerning variable $t$ to be specified, into sub-domains in order to minimize the effect of the small parameter $\varepsilon>0$. For example if we divide the domain $[1,2]$ into $10-12$ sub-domains, for $\varepsilon=0.01$ and the example addressed in the next pages, the error applying the generalized method of lines on which sub-domain is very small. Finally, we reconnect the sub-domains by solving the system of equations corresponding to the partial differential equation in question on each of these $N_{1}-1$ nodes. As $N_{1}$ is a small number for the amount of nodes (typically $N_{1}=10-12$), we have justified the terminology hyper-finite differences.

First observe that, in equation (3), $t\in[1,2]$ in $\Omega$. Discretizing in $t$ ($N_{5}$ equal pieces which will generate $N_{5}$ lines ), we recall that such a general equation (3),

[TABLE]

partially in finite differences, has the expression

[TABLE]

$\forall n\in\{1,...,N_{5}-1\}$. Here, $u_{n}(\theta)$ corresponds to the solution on the line $n$. Thus we may write

[TABLE]

where

[TABLE]

2.1 The main algorithm

To obtain a fixed point for each $T_{n}$ indicated in (13) is perfectly possible if $\varepsilon\approx\mathcal{O}(1).$ However, for the case in which $\varepsilon>0$ is small, we highlight once more the error may be relevant, so that we propose the following algorithm to deal with such a situation of small $\varepsilon>0.$

Choose $N=30-100$ and $N_{1}=10-12$ (specifically for the example in the next lines). Divide the interval domain in the variable $t$ into $N_{1}$ equal pieces (for example, for the interval $[1,2]$ through a concerning partition $\{t_{0}=1,t_{1},...,t_{N_{1}}=2\}$, where $t_{k}=1+k/N_{1}$ and $d=1/(N\;N_{1})$ is the grid thickness in $t$). 2. 2.

Through the generalized method of lines, solve the equation in question on the interval $[t_{k},t_{k+1}]$ as function of $u(t_{k})$ and $u(t_{k+1})$ and the domain shape.

To calculate $u^{k}=\{u_{n}^{k}\}$ on $[t_{k},t_{k+1}]$, proceed as follows. First observe that the equation in question stands for

[TABLE]

where

[TABLE]

Such an equation is solved through the Banach fixed point theorem, that is, defining

[TABLE]

equation (14) stands for

[TABLE]

so that for $n=1$ we have

[TABLE]

We may use the Contraction Mapping Theorem to calculate $u_{1}^{k}$ as a function of $u_{2}^{k}$ and $u(t_{k})$. The procedure would be,

1. (a)

set $x_{1}=u_{2}^{k},$ 2. (b)

obtain recursively

[TABLE] 3. (c)

and finally get

[TABLE]

Thus, we have obtained

[TABLE]

We can repeat the process for $n=2$, that is, we can solve the equation

[TABLE]

which from above stands for

[TABLE]

The procedure would be:

1. (a)

Set $x_{1}=u_{3}^{k}$, 2. (b)

calculate

[TABLE] 3. (c)

obtain

[TABLE]

We proceed in this fashion until obtaining

[TABLE]

We have obtained $u_{N-1}^{k}$ . We may then calculate

[TABLE]

[TABLE]

and so on, up to finding

[TABLE]

Thus this part of the problem is solved. 3. 3.

Calculate the solution on the lines corresponding to $u_{t_{1}},\ldots,u_{t_{N_{1}-1}}$, by solving the system,

[TABLE]

which correspond to the partial differential equation in question on the line $k\;N$, where

[TABLE]

Here may use the Banach fixed point theorem for the final calculation as well.

The problem is then solved.

2.2 A numerical example

In this section we present a numerical example. Consider the equation

[TABLE]

where, for a Ginburg-Landau type equation (see [1, 5] for the corresponding models in physics),

[TABLE]

with the boundary conditions

[TABLE]

where $\Omega=\{(r,\theta)\;:\;1\leq r\leq 2,\;0\leq\theta\leq 2\pi\}$,

[TABLE]

[TABLE]

Through the generalized method of lines, for $N=30$ (30 lines in which sub-domain), $N_{1}=10$ (10 sub-domains) and $d=1/(N\;N_{1})$, in polar coordinates and finite differences (please see [6] for general schemes in finite differences), equation (17) stands for

[TABLE]

$\forall n\in\{1,\ldots,N-1\}.$

At this point we present, through the generalized method of lines, the concerning algorithm which may be for the softwares mathematica or maple.

In this software, $x$ stands for $\theta$.

[TABLE]

[TABLE]

At this point we present the expressions for $N_{1}=10$ and $\varepsilon=0.01$. In the next lines $x$ stands for $\theta$ ($0\leq\theta\leq 2\pi$).

For each line $u[t_{k}]=U[k]$ we have obtained:

[TABLE]

Remark 2.1**.**

Observe that since $\varepsilon=0.01$ the solution is close to the constant value $1.3247$ along the domain, which is an approximate solution of equation $-u^{3}+u+1.0=0.$ Finally, the first output of the method is the solution on the $N_{1}-1=9$ nodes $U[1],\ldots,U[9]$ which, in some sense, justify the terminology hyper-finite differences, even though the solution in all the $N\cdot N_{1}=300$ lines have been obtained.

3 Conclusion

In this article we have developed two improvements concerning the generalized method of lines. For a large class of models, we have solved the problem of minimizing the error as the parameter $\varepsilon>0$ is small. In a first step we present a proximal formulation through the introduction of a parameter $K>0$ and related equation part properly specified. In a second step, we develop the hyper-differences approach which corresponds to a domain division in smaller sub-domains so that the solution on each sub-domain is obtained through the generalized method of lines.

We highlight the methods here developed may be applied to a large class of problems, including the Ginzburg-Landau system in superconductivity in the presence of a magnetic field and respective magnetic potential.

We intend to address this kind of model and others such as the Navier-Stokes system in a future research.