A Logarithmic Finite Difference Method for Troesch’s Problem

Show more

1. Introduction

Troesch’s problem arises in the investigation of the confinement of plasma column by radiation pressure and is defined by

${u}^{\u2033}=\lambda \mathrm{sinh}\left(\lambda u\right)$ (1)

subject to the boundary conditions

$u\left(0\right)=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}u\left(1\right)=1$ (2)

The main difficulty associated with Troesch’s problem is the boundary layer near $x=1$ . Accordingly, many researchers try to solve this equation numerically, some of these methods are: Sinc-Collocation method [1] , the modified homotopy perturbation technique [2] , a smart nonstandard finite difference for second order nonlinear boundary value problem [3] , the finite element method and discontinuous Galerkin methods [4] [5] , and the shifted Jacobi-Gauss collocation method [6] . A cubic spline collocation method is given in [7] . All of these methods provide a numerical solution of Troesch’s problem for moderate values of λ. Recently, an accurate asymptotic approximation of Troesch’s problem with large values up to $\lambda =60$ is reported in [8] . For more details see [9] [10] .

This paper is devoted to derive a logarithmic finite difference method as a new method to solve the Troesch’s problem for arbitrarily large parameter λ.

The paper remaining of this paper is organized as follows. In Section 2, we derive the numerical scheme as well as a fourth order finite difference scheme. Section 3 is devoted for the numerical results and comparisons with some existing methods. Concluding remarks are given in Section 4.

2. Numerical Method

In order to derive the numerical method for solving Troesch’s, we present the following definition.

Definition: The inverse sine hyperbolic function is defined by

$sin{h}^{-1}x=ln\left(x+\sqrt{{x}^{2}+1}\right)\mathrm{,}\text{\hspace{0.17em}}x\in \mathbb{R}$ (3)

and its derivative is defined by

$\frac{\text{d}}{\text{d}x}\left({\mathrm{sinh}}^{-1}x\right)=\frac{1}{\sqrt{1+{x}^{2}}}.$ (4)

By using the previous definition, the Troesch’s problem (1) can be written as

$u=\frac{1}{\lambda}{\mathrm{sinh}}^{-1}\left(\lambda {u}^{\u2033}\right)=\frac{1}{\lambda}\mathrm{ln}\left(\lambda {u}^{\u2033}+\sqrt{{\left(\lambda {u}^{\u2033}\right)}^{2}+1}\right),$ (5)

which is the main point in this paper.

Using uniform mesh including $\left(n+1\right)$ points $0={x}_{0}<{x}_{1}<{x}_{2}<\cdots <{x}_{n}=1$ , such that ${x}_{i}=ih,i=0,1,\cdots ,n$ , where $h=\frac{1}{n},{\alpha}^{-1}=\lambda {h}^{2}$ . We denote the exact and the numerical solutions respectively by $u\left({x}_{i}\right)$ and ${U}_{i}$ at the grid point ${x}_{i}$ . Now by using the second order central finite difference approximation for the second derivative

${u}^{\u2033}\left({x}_{i}\right)\simeq \frac{1}{{h}^{2}}\left({U}_{i+1}-2{U}_{i}+{U}_{i-1}\right)$ (6)

into the Equation (5), this will lead us to the logarithmic finite difference method

$\begin{array}{l}\mathrm{ln}\left[\alpha \left({U}_{i+1}-2{U}_{i}+{U}_{i-1}\right)+\sqrt{{\left[\alpha \left({U}_{i+1}-2{U}_{i}+{U}_{i-1}\right)\right]}^{2}+1}\right]-\lambda {U}_{i}=0,\\ i=1,2,\cdots ,n-1\end{array}$ (7)

subject to the boundary conditions

${U}_{0}=0\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}{U}_{n}=1.$ (8)

The resulting system in (7) represents a nonlinear tridiagonal system in the unknown solution ${\left\{{U}_{i}\right\}}_{i=1}^{n-1}$ . The details of the logarithmic can be given displayed by the following algorithm.

Algorithm 1. Logarithmic method

1) For $i=1$ using (7) and the boundary condition ${U}_{0}=0$ , we obtain

$\mathrm{ln}\left[\alpha \left({U}_{2}-2{U}_{1}\right)+\sqrt{{\left[\alpha \left({U}_{2}-2{U}_{i}\right)\right]}^{2}+1}\right]-\lambda {U}_{1}=0$ (9)

2) For $i=2,3,\cdots ,n-2$ we use

$\begin{array}{l}\mathrm{ln}\left[\alpha \left({U}_{i+1}-2{U}_{i}+{U}_{i-1}\right)+\sqrt{{\left[\alpha \left({U}_{i+1}-2{U}_{i}+{U}_{i-1}\right)\right]}^{2}+1}\right]-\lambda {U}_{i}=0,\\ i=2,\cdots ,n-2\end{array}$ (10)

3) For $i=n-1$ using (7) and the boundary condition ${U}_{n}=1$ to obtain

$\mathrm{ln}\left[\alpha \left(1-2{U}_{n-1}+{U}_{n-2}\right)+\sqrt{{\left[\alpha \left(1-2{U}_{n-1}+{U}_{n-2}\right)\right]}^{2}+1}\right]-\lambda {U}_{n-1}=0$ (11)

Newton’s method is used to solve the nonlinear system which can be described in the following way.

Algorithm 2. Newton’s method

1) For $s=0,1,2,\cdots $ .

2) Solve the tridiagonal linear system

$JZ=F\left({U}^{(s)}\right)$

for the unknown vector Z using Crout’s method. The elements of the function ${F}_{i}\left({U}_{i}\right)$ and the Jacobian matrix J are defined as follows

$\begin{array}{l}{F}_{i}\left({U}_{i}\right)=ln\left[\alpha \left({U}_{i+1}-2{U}_{i}+{U}_{i-1}\right)+\sqrt{{\left[\alpha \left({U}_{i+1}-2{U}_{i}+{U}_{i-1}\right)\right]}^{2}+1}\right]-\lambda {U}_{i}=\mathrm{0,}\\ i=\mathrm{1,}\cdots \mathrm{,}n-1\end{array}$ (12)

and the Jacobian matrix has the tridiagonal structure

$J=\left[\begin{array}{cccccc}{b}_{1}& {c}_{1}& & & & \\ {a}_{2}& {b}_{2}& {c}_{2}& & & \\ & & & & & \\ & & {a}_{i}& {b}_{i}& {c}_{i}& \\ & & & {a}_{n-2}& {b}_{n-2}& {c}_{n-2}\\ & & & & {a}_{n-1}& {b}_{n-1}\end{array}\right]$

with elements defined by

${a}_{i}=\frac{\alpha}{\sqrt{{\left[\alpha \left({U}_{i+1}-2{U}_{i}+{U}_{i-1}\right)\right]}^{2}+1}},\text{\hspace{0.17em}}i=2,3,\cdots ,n-1$

${b}_{i}=-\lambda -\frac{2\alpha}{\sqrt{{\left[\alpha \left({U}_{i+1}-2{U}_{i}+{U}_{i-1}\right)\right]}^{2}+1}},\text{\hspace{0.17em}}i=1,2,3,\cdots ,n-1$

${c}_{i}=\frac{\alpha}{\sqrt{{\left[\alpha \left({U}_{i+1}-2{U}_{i}+{U}_{i-1}\right)\right]}^{2}+1}},\text{\hspace{0.17em}}i=1,2,3,\cdots ,n-2$

3) Update your solution by using

${U}^{\left(s+1\right)}={U}^{\left(s\right)}-Z$

4) Repeat steps (1)-(3) till the following condition

${\left|{U}^{\left(s+1\right)}-{U}^{\left(s\right)}\right|}_{\infty}\le \u03f5$

is satisfied. The constant $\u03f5$ is assumed to be small.

The initial guess vector is taken as the linear interpolation between the given boundary conditions, and has the following form

${U}_{i}^{\left(0\right)}={x}_{i}=ih,i=1,2,\cdots ,n-1$

The previous method is of second order accuracy and it works for wide range values of λ.

A fourth order finite difference method can be used for solving the Troche’s method for limited values of λ. This method can be given as follows.

By using the forth order approximation of the second derivative

${u}^{\u2033}\left({x}_{i}\right)\simeq \frac{{\delta}^{2}}{1+\frac{1}{12}{\delta}^{2}}{U}_{i}\mathrm{,}$ (13)

where ${\delta}^{2}{U}_{i}={U}_{i-1}-2{U}_{i}+{U}_{i+1}$

Using this approximation and after some manipulation we will end with the following scheme

${U}_{i-1}-2{U}_{i}+{U}_{i+1}=\omega \left[\mathrm{sinh}\left(\lambda {U}_{i-1}\right)+10\mathrm{sinh}\left(\lambda {U}_{i}\right)+\mathrm{sinh}\left(\lambda {U}_{i+1}\right)\right],$ (14)

$\text{for}\text{\hspace{0.17em}}i=\mathrm{1,2,}\cdots ,n-\mathrm{1,}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{where}\text{\hspace{0.17em}}w=\frac{1}{12}\lambda {h}^{2}\mathrm{,}$ (15)

${U}_{0}=0\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}{U}_{n}=1$ (16)

The resulting method (14) is of fourth order accuracy. The numerical solution ${\left\{{U}_{i}\right\}}_{i=1}^{n-1}$ in this case can be also obtained by solving the nonlinear tridiagonal system. Newton’s method is used to solve this system as we did before.

3. Numerical Results

To test the efficiency of the proposed method, we choose $h=0.0005$ and $\u03f5={10}^{-12}$ , for the moderate values of $\lambda =0.5,1,5,10$ . Comparison with some existing methods is given in Tables 1-5. In Figure 1 we display the numerical solution for $\lambda =3,5,8,10$ . The numerical results indicate that the proposed method is highly accurate comparing to the published works.

For large values of λ, we choose
$h=0.0001$ . In Table 6 and Table 7, we display the numerical solution for
$\lambda =100,150,200,300,1000$ . To the knowledge of the authors no published work discussed these values of λ before. Again we test our method severely for
$\lambda ={10}^{4},{10}^{5}$ and 10^{6}, the method won in this case as well and the results are given in Table 8, these results are given for the first time for such values of λ up to my knowledge. We display the numerical solution for
$\lambda \mathrm{=100,200,1000,10000}$ in Figures 2-5.

Finally, we tested the fourth order method and make some comparison with the logarithmic method, we have noticed that, the results produced by this method are quite good for the small values of λ only, see Table 9, and the method is handicapped and fail for larger values of λ.

Table 1. Troesch’s problem with λ = 0.5 and h = 0.0005.

Table 2. Troesch’s problem with λ = 1.0 and h = 0.0005.

Table 3. Troesch’s problem with λ = 5.0 and h = 0.0005.

Table 4. Troesch’s problem with λ = 10 and h = 0.0005.

Table 5. Troesch’s problem with λ = 60 and h = 0.0001.

Figure 1. Numerical solution for λ = 3, 5, 8, 10.

Table 6. Numerical solution using the logarithmic method and h = 0.0001.

Table 7. Numerical solution using the logarithmic method and h = 0.0001.

Table 8. Numerical solution using Logarithmic method for large values of λ and h = 0.0001.

Table 9. Comparison of the proposed methods at λ = 0.9 with h = 0.0005.

Figure 2. Numerical solution with λ = 100.

Figure 3. Numerical solution with λ = 200.

Figure 4. Numerical solution with λ = 1000.

Figure 5. Numerical solution with λ = 10,000.

4. Conclusion

In this work, we have derived a logarithmic finite difference method to overcome the difficulty in solving the Troesch’s problem for large values of λ. This progress is very important, since all existing methods were trying to obtain the numerical solution for Troesch’s problem for large values of λ. The logarithmic method which we have derived is able and succeeds to get the numerical solution of Troesch’s problem for $\lambda =0.5,5,10,\cdots ,{10}^{6}$ . To recap things, a new logarithmic finite difference method is derived and can provide the numerical solutions for large values of λ. I think that, to the best of our knowledge, the method and some of the given results are published for the first time.

References

[1] El-Gamel, M. (2013) Numerical Solution of Troesch’s Problem by Sinc-Collocation Method. Applied Mathematics, 4, 707-712.

https://doi.org/10.4236/am.2013.44098

[2] Feng, X., Mei, L. and He, G. (2007) An Efficient Algorithm for Solving Troesch’s Problem. Applied Mathematics and Computation, 189, 500-507.

https://doi.org/10.1016/j.amc.2006.11.161

[3] Erdogan, U. and Ozis, T. (2011) A Smart Nonstandard Finite Difference Scheme for Second Order Nonlinear Boundary Value Problems. Journal of Computational Physics, 230, 6464-6474.

https://doi.org/10.1016/j.jcp.2011.04.033

[4] Deeba, E., Khuri, S.A. and Xie, S. (2000) An Algorithm for Solving Boundary Value Problems. Journal of Computational Physics, 159, 125-138.

https://doi.org/10.1006/jcph.2000.6452

[5] Temimi, H. (2012) A Discontinuous Galerkin Finite Element Method for Solving the Troesch’s Problem. Applied Mathematics and Computation, 219, 521-529.

https://doi.org/10.1016/j.amc.2012.06.037

[6] Doha, E.H., Baleanu, D., Bahrawi, A.H. and Hafez, R.M. (2014) A Jacobi Collocation Method for Tresch’s Problem in Plasma Physics. Proceedings of the Romanian Academy Series A, 15, 130-138.

[7] Khuri, S.A. and Sayfy, A. (2011) Troesch’s Problem: A B-Spline Collocation Approach. Mathematical and Computer Modelling, 54, 1907-1918.

https://doi.org/10.1016/j.mcm.2011.04.030

[8] Temimi, H. and Kurkcu, H. (2014) An Accurate Asymptotic Approximation and Precise Numerical Solution of Highly Sensitive Troesch’s Problem. Applied Mathematics and Computation, 235, 253-260.

https://doi.org/10.1016/j.amc.2014.03.022

[9] El hajaji, A., Hilal, K., El Jalila, G. and El Mhamed, M. (2014) New Numerical Scheme for Solving Troesch’s Problem. Mathematical Theory and Modeling, 4, 65-72.

[10] Scott, M (1975) On the Conversion of Boundary-Value Problems into Stable Initial Value Problems via Several Invariant Imbedding Algorithms. In: Aziz, A.K., Ed., Numerical Solution of Boundary-Value Problems for Ordinary Differential Equations.