The modified quadrature method for Laplace equation with nonlinear boundary conditions

1.   Introduction

Laplace equation arises in many important fields of science and engineering, because it is applicable to a wide range of different physical and mathematical phenomena, such as conductivity, solid mechanics, fluid, heat radiation and heat transfer. With the rapid development of computer power, many numerical methods can be used to solve the boundary value problems of Laplace equation (see Stefan et al. ^[1], Kuo et al. ^[2], Pouria et al. ^[3], Malgorzata et al. ^[4] and Li et al. ^[5]). In many numerical methods, the boundary element method may be one of the best candidates for numerical simulation due to the reduction of the dimension of boundary value problems. A lot of research on boundary integral equation methods has been devoted to developing the numerical solutions to Laplace equation with nonlinear boundary conditions. Ruotsalainen and Wendland ^[6] have used the spline Galerkin method to solve the nonlinear boundary integral equation of Laplace equation. The Wavelet-Galerkin method with Legendre wavelet functions have been used to approximate the solutions of a nonlinear boundary integral equation by Maleknejad and Mesgaran ^[7]. Pouria and Mehdi ^[8] have solved the problem numerically by the discrete collocation method using thin plate splines constructed on scattered points as basis functions. Chen, Wang and Xu ^[9] have proposed the fast Fourier-Galerkin method to solve the nonlinear boundary problem of Laplace equation. The above methods are the Galerkin method and collocation methed. The following disadvantages exist in using the Galerkin and collocation methods to solve nonlinear boundary integral equations: (1) the discrete matrix is full and each element requires the calculation of the weakly singular integral, which makes the calculation more complicated; (2) the accuracy of numerical solutions is lower ^[10,11]. However, the mechanical quadrature method was proposed by Li and Huang ^[12] to solve axisymmetric Laplace equation with nonlinear boundary conditions, which has high accuracy of the order $O(h^3)$ and low computational complexity. In this work, a modified quadrature method with the same precision and computational complexity as the mechanical quadrature method is introduced to accurately and efficiently deal with the nonlinear boundary value problem of Laplace equation, which is described as follows:

where $\Omega\subset R^2$ is a bounded, simply connected domain with a smooth boundary $\Gamma$, and $\Gamma$ is a closed curve, and the function $g(x)$ and $f(x)$ are known on $\Gamma$.

By the potential theory, (1.1) can be converted into the following boundary integral equation

where $\alpha(y)$ is related to the interior angle $\theta(y)$ of tangent lines at $y\in \Gamma$. Especially when $y$ is on a smooth part of the boundary $\Gamma$, then $\alpha = 1/2$. $k(y, x)$ is the foundation solution of Laplace equation

and

The boundary value $u$ on $\Gamma$ can be solved by (1.2), and then the normal derivative $\frac{\partial{u}}{\partial{n}}$ on $\Gamma$ can be obtained from (1.1). $u(y)$ can be calculated as follows:

Eq (1.2) is a weakly singular nonlinear boundary integral equation, whose solution exists and is unique as long as the three assumptions are satisfied ^[6]:

(1) diam$(\Omega) < 1$,

(2) $g(., u)$ is measurable for $u\in R$, and $g(x, .)$ is continuous for $x\in \Omega$,

(3) $\frac{\partial{g(x, u)}}{\partial{u}}$ is Borel measurable and satisfies $0 < l < \frac{\partial{g(x, u)}}{\partial{u}} < L < \infty$.

The modified quadrature method is presented for solving the nonlinear boundary integral equation of Laplace equation based on the composite trapezoidal rule. The modified quadrature method was first introduced by Christiansen ^[13]. Later Abou El-Seoud ^[14] proved that the numerical solutions of the order $O(h^2)$ can be obtained. Furthermore, Saranen ^[15] proved the convergence order $O(h^3)$ of the method for smooth solutions. By the modified quadrature method, the nonlinear boundary integral equation is discretized to get a nonlinear system. The calculation of the discrete matrix becomes very simple and is without any singular integrals. Moreover, the convergence of numerical solutions is proved by the theory of compact operators, and the numerical solutions of the nonlinear system is obtained by the Newton iteration. Numerical examples support the theoretical analysis.

This paper is organized as follows: in section 2, the modified quadrature method is described; in section 3, a nonlinear system can be obtained by the modified quadrature method; in section 4, according to the theory of compact operators, the convergence of numerical solutions is proved; in section 5, the Newton iteration is extensively described by the Ostrowski fixed point theorem; in section 6, numerical examples are provided to verify the theoretical results.

2.   The modified quadrature method

We consider Symm's integral equation

where $\Gamma\subset R^2$ is a closed smooth Jordan curve. Assume that $\Gamma$ be described by a $1-$periodic parameter mapping:$x(t) = (x_1(t), x_2(t)):[0, 1]\rightarrow\Gamma$ with $\mid x^{'}(t)\mid > 0$.

Defined

where the kernel $K(t, \tau) = -2\ln|x(s)-x(t)|$ and $u(t) = v(x(t))|x^{'}(t)|$. Then, (2.1) becomes

where $f(s) = g(x(s))$.

For further discussion, (2.3) is converted into the following equation

Let $h = 1/N$ be mesh width of interval $[0, 1]$ and $s_i = ih(i = 0, \cdots, N-1)$ be the nodes. Then $(Ku)(s_i)$ becomes

since the integrand $K(s_i, t)(u(t)-u(s_i))$ vanishes at the point $t = s_i$, the following approximation is obtained by the composite trapezoidal rule

For the first integral in (2.5), we get

where $x_{\rho}(s) = \rho e^{i2\pi t}$ is parametric representation of the circle with radius $\rho = e^{-1/2}$.

Since the first integral satisfies $-2\int_0^1\ln|x_{\rho}(s_i)-x_{\rho}(t)|dt = 1$ in (2.7) and the kernel of the second integral is smooth, we obtain the approximation

Replacing $u(s_i)$ with the unknown numbers $u_i$, the modified quadrature method is obtained by (2.6) and (2.8) as follows:

where $f_i = f(s_i)$.

Then, (2.9) can be converted to a matrix operator equation:

where $U = (u_0, \ldots, u_{N-1})^T$, $F = (f_0, \ldots, f_{N-1})^T$, and $B = (B_{ij})$ is $N\times N$ matrix, which satisfies^[15]

If the discretization parameter $h$ is small enough, the matrix $B$ is nonsingular, we have the error estimate^[15]

3.   Numerical procedures of the modified quadrature method

Let $\Gamma$ be described by the parameter mapping: $x(s) = (x_1(s), x_2(s)):[0, 1]\rightarrow\Gamma$, with $\mid x^{'}(s)\mid = \big[(x_{1}^{'}(s))^2+(x_{2}^{'}(s))^2\big]^{1/2} > 0$. Let $C^{2m}[0, 1]$ denotes the set of $2m$ times differentiable periodic functions with the periodic $1$ and $x_{1, 2}\in C^{2m}[0, 1]$. Define the following integral operators on $C^{2m}[0, 1]$

where $u(t) = u(x_1(t), x_2(t))|x^{'}(t)|/\pi$, $k^{*}(t, s) = k^{*}(x(t), x(s))$ is a smooth function, and $g(u(t)) = g(x(t), u(t))|x^{'}(t)|/2\pi$, $k(t, s) = -2\ln|x(t)-x(s)|$ is a logarithmic weak singular function, Then (1.2) is equivalent to

Let $h = 1/N$ be the mesh width and $t_i = s_i = ih, (i = 0, 1, \ldots, N-1)$ be the nodes.

(1) Since $K^{*}$ is a smooth integral operator, we obtain Nyström approximation with high accuracy by the trapezoidal rule.

with the error estimate

(2) Since $K$ is a logarithmic weak singular operator, we obtain Nyström approximation with high accuracy by (2.10)

where $k_h(t, s)$ is defined

By (2.12), the error estimate $K$ is

We have the numerical approximate equations of (3.3)

Obviously, (3.9) is a nonlinear equation system. When $u_h$ on the boundary $\Gamma$ are obtained, we can calculate $u_h(y), y\in \Omega$ by the following form

4.   Convergence analysis

Let's assume that the eigenvalues of operator $K^*$ and its approximate operator $K^*_h$ are not equal to 1. (3.3) and (3.9) are rewritten as

and

in which $M = (I-K^{*})^{-1}$ and $M_h = (I-K_h^{*})^{-1}$.

We use the following theorem to prove the convergence of the proposed method.

Theorem 4.1. The approximate operator sequence $\{M_h\}$ is collectively convergent to $M$ in space $C[0, 1]$, i.e

Proof. Since $k_h(s, t)$ is the continuous approximation of the integral kernel $k(s, t)$, we get ^[16] $K_h\overset{c.c}{\rightarrow}K$ in $C[0, 1]$, which shows any bounded sequence in space $\{y_m\in C^m[0, 1]\}$ has a convergent subsequence $\{K_hy_m\}$. Assume $K_hy_m\rightarrow z$ as $m\rightarrow\infty$. Because $K^{*}$ is a smooth integral operator, based on theory of compact operators ^[16], we obtain$M_h\overset{c.c}{\rightarrow}M$ in $C[0, 1]$. Further, we construct the following inequalities

where $\|.\|$ denotes the norm of $L(C^{2m}[0, 1], C^{2m}[0, 1])$. Thus we get $\{M_hK_h\}$ is compact operator sequence.

Since $K_h\overset{c.c}{\rightarrow}K$ in $C[0, 1]$, $\forall y\in C^{2m}[0, 1]$, we have $\|K_hy-Ky\|\rightarrow 0$ as $h\rightarrow 0$, and get

which shows $M_hK_h\overset{P}{\rightarrow}MK$, where $\overset{P}{\rightarrow}$ shows the point convergence. The proof is completed.

5.   Description of Newton iteration

The Newton iteration is provided to solve the nonlinear Eq (3.9). We denote

with $z = (z_0, \cdots, z_{N-1})^T = u_h$, and define

with $K^*_{ij} = K_h^*(s_i, t_j)$ and $K_{ij} = K_h(s_i, t_j)$. Then, (3.9) can be converted into the following equation

Moreover, the Jacobian matrix of $\Psi(z)$ is defined as

Thus, the New iteration can be obtained

Lemma 5.1. ^[17] (Ostrowski) Assume there is a fixed point $z^{*}\in int(D)$ of the mapping: $\omega: D\subset R^N\rightarrow R^N$ and the $F-derivation$ of $\omega$ at point $z^{*}$ exists. If the spectral radius of $\omega^{'}(z^{*})$ satisfies

Then, there is an open ball $S = S(z^{*}, \delta_{0})\subset D$ that for $z^{0}\in S$, the iterative sequence (5.5) is stable and convergent to $z^{*}$.

Lemma 5.2. ^[17] Assume $A, C\in L(R^N)$, $\|A^{-1}\| < \beta$, $\|A-C\| < \alpha$, $\alpha\beta < 1$, then $C$ is invertible and $\|C^{-1}\| < \beta/(1-\alpha\beta)$.

Theorem 5.3. Assume $\Psi: D\subset R^N\rightarrow R^N$ is $F-derivative$, and $z^{*}$ satisfies $\Psi(z) = 0$. $A: S\subset D\rightarrow L(R^N)$ is invertible and continuous at $z^{*}$, where $S$ is the neighborhood of $z^{*}$. Thus, there is a close ball $\bar{S} = \bar{S}(z^{*}, \delta)\subset S$ that $\omega$ is $F-derivative$ at $z^{*}$:

Proof. Let $\beta = \|(A(z^{*}))^{-1}\| > 0$. Because $A(z)$ is continuous and $A(z^{*})$ is invertible at $z^{*}$, when $0 < \varepsilon < (2\beta)^{-1}$, $\exists\delta > 0$, for $z\in\bar{S}(z^{*}, \delta)$, there is $\|A(z)-A(z^{*})\| < \varepsilon$. By Lemma 5.2, $A(z)$ is invertible and $\|(A(z))^{-1}\|\leq\beta/(1-\varepsilon\beta)$ for any $z\in \bar{S}$. Thus, the following function is constructed

Because $\Psi(z)$ is derivative at $z^{*}$, when $z\in\bar{S}(z^{*}, \delta)$, $\exists\delta > 0$, an inequality is obtained by the definition of the $F-derivation$ as follows:

Further, the derivation of $\omega(z)$ is considered

with $c = 2\beta(\beta\|\Psi^{'}(z^{*})\|+1)$. By the definition of the $F-derivation$, the $F-derivation$ of $\omega$ at $z^{*}$ is got.

According to the definition of $A$ in (5.4), we obtain $\rho(\omega^{'}(z^{*})) = 0 < 1$. By Lemma 5.1, the iterative sequence is stable and convergent to $z^{*}$.

6.   Numerical examples

In this section, we carry out three numerical examples for Laplace equation with nonlinear boundary conditions by the modified quadrature method, in order to verify the error in the previous sections. Let error = $|u_h-u|$, ratio = $|u_h-u|/|u_{h/2}-u|$, and $N$ denotes the mesh nodes.

Example 1. Laplace equation is considered on a plane circular domain $\Omega = \{(x_1, x_2)|x_1^2+x_2^2\leq 0.4^2\}$ with $\Gamma:(x_1, x_2) = (0.4\cos(2\pi s), 0.4\sin(2\pi s), s\in[0, 1])$. We describe the nonlinear boundary condition as: $g(x, u) = u+\sin(u)$ and $f = (x_1+x_2)+(x_1+x_2)/|x^{'}|+\sin(x_1+x_2)$. The analytic solution is $u(x) = x_1+x_2$. The numerical results are listed in Table 1.

Example 2. Laplace equation is considered on a plane elliptic domain $\Omega = \{(x_1, x_2)|x_1^2/a^2+x_2^2/b^2\leq 1\}$ with $a = 1/3$, $b = 1/2$. The boundary is described as $\Gamma:(x_1, x_2) = (a\cos(2\pi s), b\sin(2\pi s), s\in[0, 1])$. We describe the nonlinear boundary condition as: $g(x, u) = u+\sin(u)$ and $f = (a\cos(2\pi s)+b\sin(2\pi s))+(b\cos(2\pi s)+a\sin(2\pi s))/|x^{'}|+\sin(a\cos(2\pi s)+b\sin(2\pi s))$. The analytic solution is $u(x) = x_1+x_2$. The numerical results are listed in Table 2.

Example 3. Laplace equation is considered on a plane circular domain $\Omega = \{(x_1, x_2)|x_1^2+x_2^2\leq 0.4^2\}$ with $\Gamma:(x_1, x_2) = (0.4\cos(2\pi s), 0.4\sin(2\pi s), s\in[0, 1])$. We describe the nonlinear boundary condition as: $g(x, u) = u+2u^4$ and $f = 0.4\sin(2\pi s)+0.4\cos(2\pi s)+(0.4\sin(2\pi s)+0.4\cos(2\pi s))/|x^{'}|+2(0.4\sin(2\pi s)+0.4\cos(2\pi s))^4$. The analytic solution is $u(x) = x_1+x_2$. The numerical results are listed in Table 3.

From the numerical results in Tables 1–3, we can see that $\log_2ratio\approx 3$, which shows that the convergence rates of $u_h$ are at least $O(h^3)$ for the modified quadrature method.

7.   Conclusion

In this paper, the modified quadrature method is presented for the numerical solutions of Laplace equation with nonlinear boundary conditions, the innovative contributions are as follows: (1) this method was first developed to solve the nonlinear boundary integral equation with weakly kernel, and at least obtain the $O(h^3)$ order accuracy of the error, and computing entry of the discrete matrices is straightforward and simple, without any singular integrals, hence, the method is appropriate to solve weakly singular problems; (2) The convergence of this method is first proved by using the theory of compact operator, which is simpler than that of Saranen ^[15]. In the future, we plan to apply this method to solve axisymmetric problems.

Acknowledgments

The author would like to thank the anonymous referees for their valuable comments and suggestions, which led to considerable improvement of the article. This work was supported by the program of Chengdu Normal University (YJRC2018-1; CS18ZDZ02).

Conflict of interest

The author declares that he has no conflicts of interest.