Integral formulations for 1-D Biharmonic and Second Order Coupled Linear and Nonlinear Boundary Value Problems

Okey Oseloka Onyejekwe; Okey Oseloka Onyejekwe

ISSN: 2689-7636

Annals of Mathematics and Physics

Research Article Open ampcess Peer-Reviewed

Integral formulations for 1-D Biharmonic and Second Order Coupled Linear and Nonlinear Boundary Value Problems

Okey Oseloka Onyejekwe*

Author and article information

Computational Science Program, Addis Ababa University, Arat Kilo Campus, Addis Ababa Ethiopia

*Corresponding author: Okey Oseloka Onyejekwe, Computational Science Program, Addis Ababa University, Arat Kilo Campus, Addis Ababa Ethiopia, E-mail: okuzaks@yahoo.com

DOI: 10.17352/amp.000001

Received: 20 August, 2018 | Accepted: 27 October, 2018 | Published: 29 October, 2018

Keywords: Integral formulations; Boundary element method; Boundary domain interpretation; Nonlinear; Multiple differential equations

Cite this as

Onyejekwe OO (2018) Integral formulations for 1-D Biharmonic and Second Order Coupled Linear and Nonlinear Boundary Value Problems. Ann Math Phys 1(1): 001-004. DOI: 10.17352/amp.000001

Abstract

Integral formulations based on a boundary-domain interpretation of the boundary element method (BEM) are applied to develop the numerical solutions of biharmonic and second order coupled linear and nonlinear boundary value problems. The governing multiple differential equations are converted to their integral analogs by applying the Green’s identity or by double integration. The resulting integral equations are put in matrix form and solved numerically to yield both the primary dependent variable and its spatial derivative. Available benchmark solutions are applied to test the reliability of the formulation. The results are found to be in conformity with the closed form solutions and also accurately represent the physics which the problems represent.

Main article text

Introduction

Many real world problems are described by appropriate sets of differential equations which can be developed as models. These model equations are commonly differential equations which can either be linear or nonlinear. Since many of them defy analytic solutions, numerical techniques are often resorted to. Hence understanding the basic conservation laws that lead to the formulation of such problems in addition to the accompanying fluxes are required in order to offer the correct interpretation of the computed results. If however the model is found inadequate and fails to reproduce physically meaningful results or is not in consonance with the physics of the problem, then the problem formulation is revisited based on the information gathered.

Developing accurate numerical techniques that enhance this process requires an ongoing development of new ideas and techniques which can effectively provide accurate numerical solutions.in affordable computing times. From a theoretical point of view, all the conservation laws that describe physical systems lead to fluxes of the quantities conserved; for example, momentum, mass and energy fluxes. To better understand these equations, we need measurable variables such as concentration, pressure, temperature etc. This requires the use of constitutive equations which relate the fluxes to the gradient of scalar being transported. Examples of such equation are the Ficks, Fourier and Newton’s law of cooling. Once they are substituted into the governing partial differential equations, the final form of such equations are obtained in terms of measurable quantities.

Partial differential equations by their very nature deal with continuous functions and must therefore have to be discretized in space and time to arrive at numerical solutions. Discretization results in a system of ordinary differential equations (ODE), which for real world application may comprise hundreds of evolution equations that must be handled numerically. There are numerous methods for achieving this objective. Going into this is way beyond the purpose of this paper. For the purposes of the work described herein we deviate from the classical boundary element (BEM) approach by adopting two types of domain –discretized integral formulation. It is found that the resulting system of discrete equation not possess ‘local support’ as if found in finite element method (FEM) formulation but in addition possess slender coefficient matrices which are easier to handle numerically unlike the square fully populated ones that result from classical application.

Mathematical formulation

Fourth-order boundary value problems including biharmonic differential equations are of huge significance due to their several applications in such areas as applied physics and engineering. Numerical and analytical techniques have been proposed for the study of such problems (Momani and Noor) [1]. Semi-analytic techniques like the Adomian decomposition techniques including several of its variations have been applied by several authors [2,3]. A ‘non-boundary only’ approach adopted herein circumvents all numerical perplexing issues related to nonlinearity and transient scalar distribution. The integral representation of a coupled system or multiple differential equations is given as [4],

$\frac{d^{2} ϕ_{k}}{d x^{2}} = γ_{k} (ϕ_{1}, ϕ_{2}, ....... ϕ_{M}) (1)$

Where the variable k represents the equation numbering system $(1 \leq k \leq M)$ . M is he total number of equations and $γ$ is a forcing function. Equation (1) is a general representation of the 1-D Poisson equation obtained by decoupling a fourth-order differential equation. The boundary conditions at the two ends of the problem domain can be specified as:

$α_{1 k} φ_{1 k} + α_{2 k} ϕ_{1 k} = α_{3 k} k = 1, M (2 a)$

$β_{1 k} φ_{1 k} + β_{2 k} ϕ_{1 k} = β_{3 k} k = 1, M (2 b)$

A quasi -linearization of the forcing function $γ_{k}$ over each of the subdomains is initiated via the Taylor series.

$γ_{k} \approx γ_{k}^{*} + \sum_{m = 1}^{M} {{\frac{\partial γ_{k}}{\partial ϕ_{m}} |}^{*} (ϕ_{m} - ϕ_{m}^{*})} (3 a)$

The average value of dependent variable over an element at the current iteration level is $ϕ_{m}^{*}$ . Equation (1) becomes:

$\frac{d^{2} ϕ_{k}}{d x^{2}} = Λ_{1 k} + \sum_{m = 1}^{M} Λ_{2 k m} ϕ_{m} (3 b)$

Where: $Λ_{1 k} = γ_{k}^{*} - \sum_{m = 1}^{M} {{\frac{\partial γ_{k}}{\partial ϕ_{m}} |}^{*} ϕ_{m}}, Λ_{2 k} = {\frac{\partial γ_{k}}{\partial ϕ_{m}} |}^{*}$

We next introduce a weighted residual formulation for equation (3b) within a subdomain (a,b)

$\int_{a}^{b} G \frac{d^{2} ϕ_{k}}{d x^{2}} d x = \int_{a}^{b} G [Λ_{1 k} + \sum_{m = 1}^{M} Λ_{2 k m} ϕ_{m}] d x (4 a)$

We get rid of the second derivative of the dependent variable by integrating equation (4a) twice to obtain:

${[G φ_{k}]}_{a}^{b} - {[\frac{d G}{d x} ϕ_{k}]}_{a}^{b} + \int_{a}^{b} \frac{d^{2} G}{d x^{2}} ϕ_{k} d x = \int_{a}^{b} G [Λ_{1 k} + \sum_{m = 1}^{M} Λ_{2 k m} ϕ_{m}] d x (4 b)$

We define two weighting functions $G_{1} a n d G_{2}$ as $G_{1} = x - a, G_{2} = b - x$ respectively. Eq. (4a) becomes:

$G_{b} φ_{b k} - G_{a} φ_{a k} - G^{'}_{b} ϕ_{b k} + G^{'}_{a} ϕ_{a k} = \int_{a}^{b} G Λ_{1 k} d x + \int_{a}^{b} G [\sum_{m = 1}^{M} Λ_{2 k m} ϕ_{m}] d x (4 c)$

Next we approximate the dependent variable and its spatial derivative within each subdomain by introducing an osculating polynomial

$ϕ_{m} = (b - a) ξ_{1} φ_{a m} + ξ_{2} ϕ_{a m} + (b - a) ξ_{3} φ_{b m} + ξ_{4} ϕ_{b m} (5)$

Where the functions are the osculating polynomials. The integral terms in equation (4c) are evaluated to finally give the following kth equation for each subdomain.

$\sum_{m = 1}^{M} [H_{11 m} φ_{a m} + H_{12 m} ϕ_{a m} + H_{13 m} φ_{b m} + H_{14 m} ϕ_{b m}] = h^{2} Λ_{1 k} / 2 (6 a)$

$\sum_{m = 1}^{M} [H_{21 m} φ_{a m} + H_{22 m} ϕ_{a m} + H_{23 m} φ_{b m} + H_{24 m} ϕ_{b m}] = h^{2} Λ_{1 k} / 2 (6 b)$

Details of the formulation together with the element level coefficient matrices can be found in [4]. A more elegant and robust formulation of the above procedure can be initiated by applying the Green’s second identity to the stationary part of the linear diffusion operator (the so called Laplace operator) to obtain a singular integral equation. We generalize equation (1) to include convection, reaction and transient terms as well as distributed load.

$D \frac{\partial^{2} ϕ}{\partial x^{2}} - U \frac{\partial ϕ}{\partial x} - \frac{\partial ϕ}{\partial t} - μ θ = Θ (x, t) (7 a)$

D is the dispersion coefficient, $ϕ$ is the dependent variable U is the velocity in the x-direction, $Θ (x, t)$ represents an external distributed source and $μ$ is the rate constant. The auxiliary equation $d^{2} G / d x^{2} = δ (x - x_{i})$ in infinite space is used to derive the free-space Green’s function $G = | x - x_{i} | + K / 2$ with a derivative $d G / d x = 0.5 [H (x - x_{i}) - H (x_{i} - x)]$ . These together with the Greens second identity are applied to equation (7a) to obtain its integral analog

$\begin{array}{l} D [\begin{array}{l} - λ ϕ (x_{i}, t) + G^{*} (x_{2}, x_{i}) ϕ (x_{2}, t) - G^{*} (x_{1}, x_{i}) ϕ (x_{1}, t) \\ - G (x_{2}, x_{i}) φ (x_{2}, t) + G (x_{1}, x_{i}) φ (x_{1}, t) \end{array}] + \\ \int_{x_{1}}^{x_{2}} G (x, x_{i}) [U (x, t) φ (x, t) + \frac{\partial ϕ}{\partial t} + μ ϕ + Θ (x, t)] d x \end{array}$

Lagrange-type interpolations are then applied to the dependent variable and its functions within a genetic element of the problem domain. That is $ϕ \approx P_{j} ϕ_{j}$ where $P_{j}$ is the interpolation function with respect to a node j . Introducing the interpolating function into the integral equation yields the element discrete equation

$D R_{i j} + (D L_{i j} + Γ_{i k j} U_{k}) φ_{j} + T_{i j} (\frac{d ϕ_{j}}{d t} + μ ϕ_{j} + Θ_{j}) = 0 (7 c)$

To test the reliability of the above formulations, we choose a fourth order beam-deflection type differential equation whose closed form solution can easily be derived

$\frac{d^{4} ϕ}{d x^{4}} = 1 - x^{2} + e^{x} (7 d)$

$\begin{array}{l} M = \frac{d^{2} ϕ}{d x^{2}} (7 e), \frac{d^{2} M}{d x^{2}} = 1 - x^{2} + e^{x} (7 f), \\ M (0) = 0, M (2) = 4 , ϕ (0) = 0, ϕ (2) = 2 \end{array}$

$M (x) = e^{x} + 0.5 x^{2} - x^{4} / 12 - 1.5279 x - 1 , M^{'} (x) = e^{x} + x - 0.333 x - 1.5279 (7 e)$

$ϕ (x) = e^{x} + x^{4} / 24 - x^{6} / 360 - 0.254646 x^{3} - 0.5 x^{2} - 0.4204 x - 1 (7 f)$

$ϕ^{'} (x) = e^{x} + x^{3} / 6 - x^{5} / 60 - 0.8893 x^{2} - x - 0.4204 (7 g)$

The order of convergence, the $L_{2} a n d L_{\infty}$ norms of the two formulations are computed with the following formulas:

$O r d e r o f C o n v e r g e n c e = \frac{\log (E r r o r (N_{2})) / \log (E r r o r (N_{1}))}{\log (N_{1} / N_{2})} (8 a)$

where $N_{i}$ is the number of elements. By the same token, $L_{2} a n d L_{\infty}$ norms are defined as

$L_{2} = \sqrt{Δ x \sum {(ϕ_{i}^{e x a c t} - ϕ_{i}^{n u m .})}^{2}}, L_{\infty} = M a x (| (ϕ_{i}^{e x a c t} - ϕ_{i}^{n u m .}) |)$ For the purposes of comparison, we shall call the first formulation mod-1 and the latter mod-2. The numerical results are hardly distinguishable for this particular example. However Table 1 shows that mod-2 displays slightly better accuracy and convergence.

We put the robustness of mod-2 to test by solving a beam deflection problem in structural mechanics. A supported beam carries a uniform load of intensity $w_{0}$ and a tensile force N. The governing fourth order differential equation for the system is given by:

$\frac{d^{4} ϕ}{d x^{4}} - \frac{N}{E I} \frac{d^{2} ϕ}{d x^{4}} = \frac{w_{0}}{E I} (8 b)$

Where EI is the bending rigidity and $ϕ$ is the displacement. The boundary conditions are

$ϕ = d^{2} ϕ / d x^{2} = 0 a t x = 0 a n d x = L (8 c)$

Changing the variables to $υ = x / L, y = E I ϕ / w_{0} L^{4}$ converts equation (10a) into

$\frac{d^{4} ϕ}{d υ^{4}} - ς \frac{d^{2} ϕ}{d υ^{2}} = 1, ς = \frac{N L^{2}}{E I} (8 d)$

With the following boundary conditions:

${ϕ |}_{υ = 0} = {\frac{d^{2} ϕ}{d υ^{2}} |}_{υ = 0} = {ϕ |}_{υ = 1} = {\frac{d^{2} ϕ}{d υ^{2}} |}_{υ = 1} (8 e)$

We adopt the same decoupling procedure in example (1). Table 2 shows the displacement profiles and their gradients across the beam. Numerical results justify the relative signs of the variable $ς$ . Identical values are recorded at the middle of the beam.

Our next example deals with transient nonlinear fourth-order equation that is described by the following equation:

$\frac{\partial ϕ}{\partial t} + χ \frac{\partial^{4} ϕ}{\partial x^{4}} - \frac{\partial^{2} ϕ}{\partial x^{2}} + Θ (ϕ) = 0, x \in [a, b], t \in [0, T] (9 a)$

Where $Θ (ϕ) = ϕ^{3} - ϕ$ , with initial and boundary conditions:

$ϕ (x, 0) = ϕ_{0} = x^{2} {(1 - x)}^{2}, ϕ (0, t) = ϕ (1, t) = 0, ϕ_{x x} (0, t) = ϕ_{x x} (1, t) = 0 (9 b)$

Equation (11) is decoupled according to:

$M (x, t) = ϕ_{x x} (x, t), M (0, t) = M (1, t) = 0 (9 c)$

$\frac{\partial ϕ}{\partial t} + γ \frac{\partial^{2} M}{\partial x^{2}} - M + Θ (ϕ) = 0, ϕ (x, 0) = ϕ_{0} = x^{2} {(1 - x)}^{2}, ϕ_{x x} (0, t) = ϕ_{x x} (1, t) = 0 (9 d)$

Equation (9) is solved with mod-2 and the nonlinearity is resolved by the Picard algorithm. Figures 1,2, show how the nonlinear initial profile decays with time

Conclusion

Two integral formulations based on a boundary-domain interpretation of the boundary element method have been used to numerically solve a Biharmonic and second order coupled linear and nonlinear boundary value problems. Unlike the BEM both methods require that the domain of the problem be discretized in the space-time domain. As a result of this domain- localized integral approach, local support is obtained among the computational nodes, and the resulting coefficient matrix is slender and amenable to numerical manipulation.

References

Momani S, Noor MA (2007) Appld Math and Computation 191: 218.
Wazwaz AM (2002) Jnl Comput Math 79: 345.
Momani S, Moadi K (2006) Appld Math and Comput 22: 185.
Ramachandran PA (1994) Boundary Element Methods in Transport Phenomena, Comput Mech Publ Elsevier Appld Sci. London, N.Y.

Copyright

© 2018 Onyejekwe OO. This is an open-ampcess article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

View Similar Articles

Share your thoughts and experiences

Order for reprints

Article Alerts

Subscribe to our articles alerts and stay tuned.

This work is licensed under a Creative Commons Attribution 4.0 International License.

Quick Enquiry

Table 1: Error Norms and Order of Convergence for mod-1 and mod-2.
No. Elements	norm (mod-1)	norm (mod-2)	Order of Convergence (mod-1)	Order of Convergence (mod-2)	norm (mod-1)	norm (mod-1)
10	.012254	.008652	1.89752	1.99658	0.01287	.00231
20	.011536	.006592	1.97321	2.00791	0.00322	.00182
40	.0083712	.001229	2.01376	2.00362	.001462	.00108

Table 2:
Coordinates	Displacement	Gradient	Displacement	Gradient
0.1	-0.5487e-01	-0.4467e+00	-0.3269e-01	-0.3075e+00
0.2	-0.9072e-01	-0.2754e+00	-0.6155e-01	-0.2666e+00
0.4	-0.1185e+00	-0.4476e+00	-0.1042e+00	-0.1518e+00
0.5	-0.1156e+00	0.7466e-01	-0.1156e+00	-0.7466e-01
0.65	-0.9573e-01	0.1850e+00	-0.1162e+00	-0.7466e-01
0.8	-0.6155e-01	0.2666e+00	-0.9072e-01	0.2754e+00

Annals of Mathematics and Physics

Integral formulations for 1-D Biharmonic and Second Order Coupled Linear and Nonlinear Boundary Value Problems

Okey Oseloka Onyejekwe*

Author and article information

Abstract

Main article text

Introduction

Mathematical formulation

Conclusion

References

Copyright

View Similar Articles

Share your thoughts and experiences

Article Alerts

Table of Contents

Submit your next article Peertechz Publications, also join of our fulfilled creators. Submit a Manuscript

© Peertechz Publications Inc., 10880 Wilshire Blvd., Suite 1101, Los Angeles, California, 90024, USA