On the Numerical Solution of Boundary Value Problem (BVP) of the Ordinary Differential Equation (ODE) - The Case of Steady-State Bio-Heat Equation with Combined Heat Transfer Coefficient by Pseudo-Spectral Collocation Method

—Spectral methods for the solution of a boundary value problem of an ordinary differential equation are reviewed with particular emphasis laid on pseudo-spectral collocation method. The pseudo-collocation method is then used to solve the one dimensional bio-heat equation with metabolic heat generation in cylindrical coordinates applied to human tissue. It was noticed that an increase in heat transfer coefficient (hA), enhanced the temperature but a decrease in the tissue thickness was observed when this coefficient was increased. The effects of the combined heat transfer coefficient are analyzed and the results indicate that the obtained solution can be used in the study of the thermal behaviour of a biological system with the potential to locate tumours in the living tissue.

Abstract-Spectral methods for the solution of a boundary value problem of an ordinary differential equation are reviewed with particular emphasis laid on pseudo-spectral collocation method.The pseudo-collocation method is then used to solve the one dimensional bio-heat equation with metabolic heat generation in cylindrical coordinates applied to human tissue.It was noticed that an increase in heat transfer coefficient (hA), enhanced the temperature but a decrease in the tissue thickness was observed when this coefficient was increased.The effects of the combined heat transfer coefficient are analyzed and the results indicate that the obtained solution can be used in the study of the thermal behaviour of a biological system with the potential to locate tumours in the living tissue.

I. INTRODUCTION
T HIS article is organized as follows; section one is the introduction to the concept of mathematical modelling in many aspects of human life, the various methods that are used to find the solutions of the models, their merits, and demerits.A brief literature review on the solution of the bioheat equation is presented.In section II, spectral methods as applied to one class of differential equations-boundary value problems (BVPs) of ordinary differential equations (ODEs) are reviewed.Chebyshev polynomials and some of their properties are introduced.Section III makes application of the pseudo-spectral method(collocation) to the solution of the steady state bio-heat equation with metabolic heat generation.Section four is the discussion of results and section five is the conclusion.
Research work done in the fields of engineering, physics, medicine, biology, and economics among others requires in some way mathematical modelling of some physical phenomena i.e., real world issues represented in mathematical terms.For us to understand and predict the behaviour of some aspects of these phenomena, a 'picture' of the physical system with simplifying assumptions is made into a mathematical framework that we can then seek to find a solution.The dynamical aspects of these phenomena are usually described by a set of differential equations both ordinary differential equations (ODEs) and partial differential equations (PDEs) giving rise to a mathematical model.As a simple example in physics, the temperature of a cup of tea in a room of constant temperature will cool over time at a rate proportional to the difference between the room temperature and the temperature of the tea.We can model this by Newton's law of cooling which is a first order ordinary differential equation, we obtain its solution and use it to predict the temperature of the tea at any given time.
There is a wide range of methods both analytical and numerical that have been developed over time for the solutions of both ODEs and PDEs.Analytical solutions are limited in their use, especially for non-linear differential equations and high order equations.For this reason, numerical solutions are preferred.Techniques for obtaining numerical solutions of differential equations are numerous in the literature, for example, see [ref].Finite difference method (FDM), finite element method (FEM), the shooting method, and spectral methods among others are some of the numerical methods commonly employed.FEM uses lower order interpolants usually linear interpolation.As a comparison, unlike in FDM, in spectral methods, the value of a derivative at a certain point in the domain space depends on the solution at all other points in space and not just on neighbouring grid points.Due to this aspect of the spectral methods, they have very high order of approximation.This is called spectral convergence in which the error decreases exponentially with an increase in grid resolution.Therefore for spectral methods, the rate of convergence depends only on the smoothness of the solution.The concept of smooth functions mean for example, C ∞ (D) : class of functions f : D to R having infinitely many continuous derivatives.Another advantage that spectral methods have over FDM is that, once the approximate spectral coefficients have been found, the approximate solution can be evaluated at any point in the domain, whereas to evaluate a finite difference solution at an intermediate point requires a further step of interpolation.Both FEM and FDM methods handle very well two dimensional space problems.
Spectral methods approximate the solutions of differential equations by means of a truncated series of orthogonal functions.Spectral methods derive their strength from a class of discretization methods called the Method of Weighted Residual (MWR) which uses the trial (or basis) functions and the test functions that ensure the approximate solution of the differential equation is as close to the exact solution as possible.Spectral methods can be applied to the following classes of differential equations.Boundary Value Problems (BVPs) for both ODEs and PDEs, Initial Boundary Value Problems (IBVPs) for PDEs (parabolic and hyperbolic), both linear and non-linear problems, the group of linear BVPs include singularly perturbed problems and eigenvalue problems.
There are three types of spectral methods, pseudocollocation, Tau and Galerkin.The choice of the method depends on the application.Collocation method deals mainly with non-linear problems/complicated coefficients.The Tau method handles complicated boundary conditions(where Galerkin is impossible).Galerkin method is used when optimal error estimates are required and also provides convenient analysis.The standard approach for application of spectral methods is to obtain the trial and test functions.There are disadvantages however, (i) the matrices resulting from discretization process have increased condition number meaning that rounding errors reduce the expected theoretical exponential accuracy, (ii) discretization matrices are not sparse so efficient algebraic solvers are difficult to apply.The above disadvantages become obvious when solving 4th order problems where stability and numerical accuracy are lost when applying higher order approximations.One solution to these challenges is to address the flexibility in the choice of test and trial functions, and reduce condition number and band-width of matrices.
The heat transfer in a living biological tissue involves a combination of three major processes of heat transfer, namely, thermal conduction in tissues, convection and perfusion of blood, and metabolic heat generation.This was first noted by Harry Pennes, an American physician and clinical researcher in 1948.The bio-heat equation emerged from this observation [1].It is also called Pennes bio-heat equation, it relates the rate of heat transfer between the blood and tissue as being proportional to the product of the volumetric perfusion rate(the volume flow rate of blood to the tissue per unit volume tissue) and the difference between the arterial blood temperature and the local tissue temperature.This equation is the typical model used for predicting the distribution of temperature inside a biological system.Pennes validated his model against a series of experimental studies he performed.It turned out that the results of Pennes bio-heat model were in reasonable agreement with the experimental data.Since then both analytic and numerical solutions of the bio-heat equation have been reported by numerous researches.Some of them are reviewed here.Zhou et al. [4] presented a steady-state one dimensional solution in Cartesian coordinates.A similar analysis in cylindrical geometry showing the parameter effects in spatial temperature distribution had been done earlier by [4].Shih et al. [5] discussed the analytic solution of living tissue model with sinusoidal heat flux on the skin surface.The obtained solution can be used to analyse the of effects metabolic heat generation, tissue thermal conductivity, blood perfusion, and heat transfer coefficient.Heat transfer within a perfused tissue in the presence of a blood vessel was modelled by Huang et.al [8].Exact analytical solution of the bio-heat equation subjected to intensive moving heat source was examined by Mohammad et.al [9].They used eigen function expansion to find the solution with potential to be used as a verification tool for numerical solutions.Hongyun et.al [10] obtained an analytical solution of the one dimensional Pennes bio-heat equation by Fourier transform and applied it to the case of multiple electromagnetic heating pulses.In the present study, we obtain the solution of the one dimensional bio-heat equation in a cylindrical tissue using the pseudocollocation method.

II. NUMERICAL SOLUTION OF AN ORDINARY DIFFERENTIAL EQUATION
Consider a boundary value problem(BVP) of an ordinary differential equation (ODE) of the form: where L is a differential operator(linear or non-linear), u = u(x) is the solution function we are looking for over the domain [−1, 1].f is a known source function( may depend on x or is a constant), u(−1) = 0 and u(1) = 0 are the boundary conditions.We denote the approximate solution function to the BVP by u N (x) and write it as an expanded finite series of orthogonal polynomials, where N + 1 is the finite number of basis functions that are considered.The commonly used orthogonal polynomials include Chebyshev, Legendre, and Bessel.Hence we write where {a k } N k=0 are the expansion coefficients we are going to determine and {ϕ k } N k=0 are basis or trial functions(orthogonal polynomials).Form the quantity called the residual.Since the basis functions are given, the challenge is to keep R[u N (x)] as small as possible across the domain [−1, 1] as we find the expansion or spectral coefficients {a k } N k=0 .
A. Determination of Expansion Coefficients {a k } N k=0 : There are three popular methods that can be used to find these coefficients.Here we describe them but use pseudo-spectral collocation method in the solution of the bioheat equation in our application.
(i) Tau-Lanczos method: {a k } N k=0 are selected such that the boundary conditions are satisfied identically and R[u N (x)] is orthogonal to as many basis functions as possible (ii) Galerkin method: The first basis functions are recombined {ϕ k (x)} N k=0 → {Φ k (x)} N k=0 so that the boundary conditions are satisfied identically.Then, the coefficients {a k } N k=0 are determined such that the residual R[u N (x)] be orthogonal to as many new basis functions {Φ k (x)} N k=0 as possible (iii) Collocation method: {a k } N k=0 are selected such that the boundary conditions are satisfied.The rest of the coefficients are determined so that R[u N (x)] vanishes at as many spatial locations as possible.The spatial locations are called collocation points.The best choice of collocation points are the Gauss-Lobatto points x k = − cos( πk N ), k = 0, 1, 2, . . ., N

B. Chebyshev Series
Fourier series is a good choice for periodic functions.Problems with non-periodic boundary conditions require trial functions based on orthogonal polynomials.A popular choice is Chebyshev polynomials defined on −1 ≤ x ≤ 1.They are defined as follows: from trigonometry cos(0θ) = 1, cos(1θ) = cos θ, cos(2θ) = 2 cos 2 θ − 1, cos(3θ) = 4 cos 3 θ − 3 cos θ, cos(4θ) = 8 cos 4 θ − 8 cos 2 θ + 1, . . .Hence cos(nθ) can be expressed as a polynomial in cos θ i.e cos(nθ 2 ) so cos(nθ) + cos(n − 2)θ = 2 cos θ cos(n − 1)θ from which we deduce that subsequent polynomials can be found by using the following recurrence relation A function u(x) is approximated via a finite series of Chebyshev polynomials as where a k are the N + 1 Chebyshev coefficients.A function interpolated by higher order polynomials leads to oscillations at the end points of the domain(Runge phenomenon) [6].If a function is interpolated on an equidistant grid, the error grows by 2 N .However, using non-equidistant distribution of points with denser points towards the domain boundaries, it can be shown that interpolation errors decrease exponentially.A common distribution of points for Chebyshev polynomials are the Gauss-Lobatto points.: C. Some Properties of Chebyshev Polynomials • boundary conditions t k (1) = 1, and t k (−1) = (−1) k

III. SOLUTION OF THE BIO-HEAT EQUATION
Consider the second order BVP of an ODE in the form of steady-state Bio-Heat Equation in cylindrical coordinates with metabolic heat generation on the domain D = [0, R] where R is the radius of the cylindrical tissue.
where Q m is the metabolic heat generation in the tissueassumed to be homogeneously distributed throughout the tissue, T a is the temperature of the arterial blood, T is tissue temperature, k is the thermal conductivity of the tissue assumed to be uniform, r is the radial position in the limb, M blood perfusion rate parameter.We use a convective boundary condition at the skin surface where h A is the combined convection/radiation heat transfer coefficient between the skin surface and the surroundings, which are at temperature T ∞ .The approximate numerical solution for the problem is cast as follows where t * k (x) is the shifted Chebyshev polynomial of the first kind and {a k } N k=0 are the spectral coefficients to be determined.Note that the interval [0, R] transforms as follows: . Equation (1) may be written as where L = kr d 2 dr 2 + k d dr − M r is the differential operator.For the bio-heat equation, let us take N = 4 which will yield a system of (N + 1) = 5 linear equations.i.e The residual is given by A. Collocation(pseudo-spectral) method: boundary conditions in Equation ( 7), we have The other three relations can be evaluated in order to solve for all the coefficients.The Chebyshev grid points on the domain [−1, 1] are given by r j = − cos( πj 4 ), j = 0, 1, 2, 3, 4, or We now work out solution for the three interior points.It can be shown from Equation (4) that From which we can find the first derivative with respect to r dT 4 (r) and the second derivative with respect to r For j = 1 and therefore , Equation ( 6) may be written as yielding where the c ′ s are the coefficients of the resulting linear system of equations for the spectral coefficients a 0 , a 1 , a 2 , a 3 and a 4 .
Similarly, for j = 2, r 2 = R 2 and j = 3, with We now solve the system of linear equations, (8), ( 9), ( 14), ( 15) and (16) simultaneously using mathematical software.A mathematical algorithm for the c ′ s is incorporated in the a code in order to enable variation of parameters during computer simulation.

IV. NUMERICAL RESULTS AND DISCUSSION
This section discusses the effects of the combined heat transfer coefficient in a living tissue using pseudo-spectral collocation method solution.Note that we used N = 4 which is equivalent to taking the first N + 1 = 5 basis functions(Checbyshev polynomials) in the solution.To determine the effects of combined heat transfer coefficients, various values are taken for observation.Figure 1 above, represents the graph of tissue thickness against body temperature using pseudo-spectral collocation method.

V. CONCLUSION
In many practical applications, the heat transfer rate in a living tissue is vital since it influences the general body temperature.The present work, helps us in understanding numerically as well as physically the effects of combined heat transfer coefficient in a living tissue.The results indicate that the obtained solution can be used in the study of the thermal behaviour of a biological system.Thus, applications of the effects of the combined heat transfer coefficient are recommended for locating tumours in living tissue in medicine.Based on the results, the effects of increasing the combined heat transfer coefficient in a living tissue which had significant effect on temperature and tissue thickness were as follows; From the numerical results, the positive values of heat transfer coefficient, hA > 0 are utilised in our computations.This corresponds to the heat transfer problem with respect to the application.The heat transfer problem is often encountered in various biological system applications for example in locating tumours in the living tissues.
It was noticed that an increase in heat transfer coefficient (hA), enhanced the temperature but a decrease in the tissue thickness was observed when this coefficient was increased.
On the Numerical Solution of Boundary Value Problem (BVP) of the Ordinary Differential Equation (ODE) -The Case of Steady-State Bio-Heat Equation with Combined Heat Transfer Coefficient by Pseudo-Spectral Collocation Method Benard A. Odongo, Alfred W. Manyonge, Dancun O. Owego, Richard O. Opiyo, and Thomas M. Onyango

Fig. 1 .
Fig. 1.The Graph of tissue thickness against body temperature.

TABLE I :
Thermo-Physical Properties of the Human Body Tissue