New Parallel Algorithms for BVPs in ODEs

صخلملا  لئاسـم لـحل ةيزاوتم ةديدج ةقيرط ريوطت وه يسيئرلا ثحبلا اذه فده لضافتلا تلاداعملا يف ةيموختلا ميقلا ةيزاوتم تابساح يف ذيفنتلل ةمئلام ةيدايتعلاا ةي عون نم MIMD ) دحاو نآ يف ةددعتم تايلمع تاذ تابساح .( ةيـصاخ ةسارد تمت ةداـيز قـيرط نع ةقيرطلا نيسحت متو ةقيرطلل أطخلا ىلع ةرطيسلاو ةيرارقتسلأا يعجرلا لماكتلاو تانيمختلا ددع . ةـفاجلا ةـيموختلا ميـقلا لئاسم ةجلاعم تمت امك ةركتبملا ةقيرطلاب .  ABSTRACT  The main objective of this paper is the development of a new parallel integration algorithm for Solving Boundary Value Problem (BVPs) in Ordinary Differential Equation, (ODEs). This algorithm is suitable for running on MIMD computing systems. We will analysis the stability and error control of the developed algorithm .We shall also consider the treatment of stiff boundary value problems by developed technique


Introduction
Numerical solution of BVPs in ODEs is a very active research area. Mostly the numerical solution of these problems takes a lot of computer time, specifically, if the integration interval is very large or if the system describing the problem consists of a large number of differential equations or if the problem is a stiff problem (consult the following references for more detail Cash (1995), Khalaf (1988Khalaf ( , 1990, Khalaf and Al-Wajih (2000), Khalaf and Al-Murshid (2000). Hence the objective of this research is the development of a new parallel algorithm which combines parallel integration processes with parallel interpolation to estimate the unknown BCs.
MIMD Computer Consists of several processors, each processor has its own memory and processing unit. These processors communicate through a suitable communication network. (for more detail, see Meiko (1989), Khalaf (1990), Al-Wajih (1999) and Al-Murshid (2000)). In MIMD computer each processor can carry out its own set of instructions, often on its own set of data, independently of all the other processors. Such computers usually number their (more complex) processors in tens rather than thousands that may be found in SIMD (Single Instruction Stream with Multiple Data Stream) computers. MIMD computers are well suited to algorithmic parallelism in which problems can be separated into concurrent independent processors (Brocklehurs, 1992).

1-The Numerical Methods for BVPs:
The well -known numerical methods for solving BVPs in ODEs are: i) Finite difference methods (Fd), ii) shooting methods and iii) collocation methods. Finite difference methods based on dividing the given interval of the independent variable by nodes and then approximating the differential equation by a given finite difference formulas at each node and this will produce a set of algebraic equations mostly non -linear which can be solved by Newton iteration or one of its alternatives (for more detail see Khalaf (1988)). To get accurate results by using these methods we have to increase the number of nodes which will produce a great number of algebraic equations which increase the complexity of the solution and takes a lot of computer time. Mostly, the iteration process will diverge because the approximation processes at the nodes will create noisy data and this noise can be accumulated by iteration processes and render the solution meaningless.
Shooting methods based on dividing the integration interval to n subintervals. At the beginning of each subinterval, values estimated for the given dependent variables then the ODEs of the problem integrated in the subinterval, by using the estimated values and then at the end of each subinterval the corresponding estimated and integrated values of corresponding depended variables are matched, i.e. at the end of each subinterval matching functions are defined. The estimated values can readjusted by Newton iteration or one of its alternatives. The problem with these methods is that the estimated values need to be very close to real solution, otherwise the iteration processes will diverge. (For more detail see for example Keller (1988), Khalaf (1988Khalaf ( , 1990). For Collocation methods see Asher et al. (1.988), and Khalaf (1997)) based on approximating the solution of the ODE by a linear combination of a set of independent simple functions. The coefficients of the combination can be estimated ( using the boundary conditions BCs and substituting the differentiation of the approximate solution at given nodes of the independent interval) by Newton methods or by some other iterative methods. Hence the iteration will diverge at most if there are noises or error in combination coefficients so that we try here to develop a new method for solving BVDs in ODEs which uses integration's and interpolation instead of iteration processes. The new methods are more suitable for running on MIMD computers, (see KLhalaf and Matti (2001) and Khalaf and Sawoor (2002)

2-The new method:
Let us consider a two point BVP (this will not affect the generality of the method): The above problem can be reduced to the following system.
To integrate system (2.2) in the interval (a,b) we need a value of y 2 (a) which is unknown. To get this value, we proceed as follows: Process 1: estimate a value S 1 for y 2 (a) and integrate the system (2.2) in the interval (a,b), we get y 1 (b) =-m 1 . Process 2: estimate another different value 82 for y2(a) and integrate the system (2.2) in the interval (a,b), we get y 1 (b) =m 2 , and so on Process n: estimate another different value Sn for y 2 (a) and integrate the system (2) in the interval (a,b) , we get y 1 (b) =m n . Hence we get the following table of data: we can write the above table as: Then by using Newton divided difference formula for (a) corresponding to 2 tion,we can find the value of y interpola And the approximate value of y 1 (a) corresponding to y 1 (b)= β is y 2 (a)=p(β)

3-ParalleIising of the Algorithm
It is clear from the last section that the integration processes 1,2,3,...,n are independent processes. Hence each integration process can be carried out in a processor, i.e. each integration process can be assigned to one of the processors of a MIMD computer, which consists of p processors. If p>n then integration processes can be speeded up by S p =n, where S p is the speed-up factor and defined as follows: if p<n and n/p=r (r is a positive integer) then the integration processes can be speeded-up approximately by the factor S p =r. It is also clear that finding Δ.S o , Δ.S 1 , ...,Δ.S n-1 are independent processors, so that they can be calculated simultaneously. Calculation of Δ 2 S o , Δ 2 S 1, … Δ 2 S n-2 are independent so that they can be calculated in parallel and the higher order calculation Δ k S i , i=0,l,...,n-k, k=3,...,n, S i =S(m i ) can be done similarly.

4-Testing the new method
In this section we try to show the effectivenesses of the new method for solving linear BVPs and non -linear BVPs: To integrate the above system by the new method, we give the following values for y 2 (0): -1,0,2,3, and we integrate the system by Euler method using step size h=0.5 we get the following table: We rewrite the above table as: Where the exact value of y2(0) is 1.
The accuracy of the final value of y 2 (0) can be increased by reducing the step size and increasing the number of estimations or the accuracy can be increased by using higher order integration methods. Example (2): (a non -linear problem) Consider the following nonlinear BVPs: y"=6 3 1 y , y(0)=0, y(l)=l, y(x) = x 3 To integrate the system by the new method we need estimation for the unknown BC which is y'(0). The result obtained by the method for y'(0) is 0.671 and the exact value is 0. The accuracy of the result can be increased by: i) reducing the step size h. ii) by increasing the number of estimation and iii) by using a higher order integration method.

5-Improving the method by using higher order integration methods
To solve example (1) by the method we give the following value for y2(0): -1,0,2,3 using forward integrating process but we integrate the system by a higher-order method, such as, fourth order R-K method using step size h=0.5 ,we get the following results: We rewrite the above tables as: And since the exact value of y 2 (0) is 1, i.e. the error of this method is 4.29817×l0 -2 . So that the effectiveness of a higher-order integration method is clear.

6-Improving the method by backward integration techniques:
To improve the performance of the parallel integration and interpolation method for BVPs , we have used Backward integrating processes to estimate the value of unknown boundary condition (BC) y2(a) . This improvement is done by integrating the system from b to a by using negative step size -h , h > 0 . The 2nd -order boundary -value problem (BVP): y"=f(x,y,y') . a  x  b (6.1) β , y(b)=  y(a)= Can be reduced to a system of equations of the form : =f(x,y To integrate system (6.2) in the interval (a,b) we need a value of y 2 (a) Which is unknown. To get this value we proceed as follows: Process (I): estimate a value S 1 for y 2 (b) and integrate the system (6.2) from b to a by using negative step size-h , h > 0 i.e. using one of the single-step method such as Euler method, we get: y 2 (a) = r 1 Process (2): estimate another different value S 2 for y 2 (b) and integrate the system (6.2) from b to a by using negative stepsize -h,h>0 i.e. using one of the single-step method such as Euler method, we get y 2 (a) = r 2 ,and so on Process (n): estimated another different value Sn for y 2 (b) and integrate the system (6.2) from b to a by using negative step size -h,h>0 i.e. using one of the single-step method such as Euler method, we get y 2 (a) = r n .
Hence we get the following table of the data: Then by using Newton divided difference formula for interpolation we can find the value of y 2 (a) corresponding to y 1 (b) = β. Here we do not need the conversion of the table.

6-1 Parallelising of the Backward integration algorithm
It is clear from the last section that the integration processes 1,2,3,..., n are independent processes. Hence each integration process can be carried out in a processor, i.e. each integration process can be assigned to one of the processor of a MIMD computer that usually number their processors in tens rather than thousands that may be found in single instruction stream with multiple data stream computers, which consists of p processors, if p>n then integration processes can be speeded-up by Sp=n, where Sp is the speed-up factor and defined as follows: S p = If p<n and n/p=m (m is a positive integer) then the integration processes can be sped up by the factor S p =m. It is also clear that finding Δ r o , Δ r 1 ,..., Δ r n-1 are independent processors, so that they can be calculated simultaneously. Calculations of Δ 2 r o , Δ 2 r 1 ,..., Δ 2 r n-1 are independent so that they can be calculated in parallel and the higher order calculations Δ k r i , k=3,...,n, r i = r(S i ) can be done similarly.

Testing of the newly developed method Example(4):
Let us consider the BVP of example (1). y"=-y'+2y y(0)=l, y(l)=e, y(x)=e x which can be reduced to the following systems 1 y =y2 y1(0)=1 2 y =-y2+2y, y1(1)=e To integrate the above system by using backward integrating process, we give the following values for y 2 (b)=y 2 (l)= -1,0,2,3 and integrate the system by Euler method using negative step size -h,h=0. The value of y2(a) by using forward integrating process is 0.7548172 and the exact value is 0. To integrate the system by the newly developed method that uses backward integrating process: first, reduce the BVP to the following system: x 1 =x o -h= 1.5 x 2 =a=x 1 -h= 1.0 By estimating the following values for y 2 (b)=y 2 (2)= -1,0,2,3; and integrating the system by Euler method with negative step size -h, h=0.5>0, we get the following results:  P(β)=y 2 (a)-y 2 (l)=0.1597225, this value which is obtained by backward integrating process is more accurate than that obtained by forward integrating process , y 2 (a)= 0.7548172, compared with the exact value 0.

6.3-Improving the performance of the backward integration process:
The accuracy of the results that are obtained by backward integration process can be increased by increasing the number of estimations provided that Δ 2 r i , Δ 3 r i , …, Δ k r i , k=2,... ,n-l are not equal to zero.

7.Stability and Convergence Analysis of the Method:
It is clear that the developed method consists of two stages: stage of integration and stage of interpolation. Stability of the integration stage is well discussed by many authors and researchers such as Lambert (1974), Henrici (1962), Gear (1971) Conte and de Boor (1981) and others.
To control the induced instability, we have to use methods of integration whose finite difference equations are of the same order as the order of differential equation (for more detail see for example Khalaf (2000). To control the inherent instability we have to use multiple shooting or by using backward integration (for more detail see khalaf (2000)).
Another way for controlling the phenomena of instability of the numerical integration methods is improving the stability interval of the methods by using the principle of the fixed point iterations. This technique is used successfully by Abdulah (2000, 2001)) and then used by Khalaf and Mahmood (2001) for developing new methods for solving chaotic systems. It is well known that (see Lambert (1974)) Stability + consistent convergence of the integration method, hence we can control the convergence character of the integration stage very easily.
We have no stability problem in the stage of interpolation. Hence there is no convergence problem in this stage. We only need the interpolation to be as accurate as possible by increasing the number of estimations and using higher order polynomials to interpolate the unknown boundary condition, i.e. we need in interpolation stage to reduce the error of interpolation. This can be done easily because the error of the interpolation at tlie point x is given by Conte and de Boor (1981):   (a,b). It is clear from the form of e n ( x ) that the error of interpolation decreases as the number of estimations increases. Hence, we can increase the accuracy of the method by increasing the number of estimations and using accurate convergent integration method. Now we can state the following: The new method is convergent if the integration Proposition: stage is convergent and the interpolation stage is accurate.

8-Treatment of Stiff Boundary Value Problem by the New Methods:
Let us consider the following Boundary value Problem: Use implicit methods such as trapezoidal method, whose region of stability is the entire negative half-plane, so that h is unrestricted by the stability requirement (see Gear (1971)). But solving the problem by trapezoidal method will lead to a system of algebraic equations which can be solved by fixed-point iterations or by Newtonian iteration, which increases the complexity of the method as well as these two techniques converge under hard conditions. So to be out of the new problem we can modify the explicit Runge-Kutta method of order 4, which that we can use a larger h. The standard Runge-Kutta of order 4 is: The new iteration formula of Runge-Kutta becomes: Using this new formula for solving the above stiff problem we require h to be 1000h<2.8(l+a). In case a=0.9, h<0.00532, which is about two times larger than that used by standard Ruuge-Kutta. Any how using oc»l, will change the convergence requirement of the method (for more detail the reader can consult the following references Khalaf and Abudullah (2000), (2001)).

9-Conclusions:
We have developed a new parallel integration algorithm for solving boundary value problems suitable for MIMD computing systems. We showed how the results of the method can be improved by increasing the number of estimations and using backward parallel integrations.
The stability and error control of the method are well analysed. The treatment of the stiff boundary value problems is considered and we have developed explicit integration algorithms for solving these stiff boundary value problems.