Numerical solution of ordinary differential equations pdf

New exact solutions were obtained for several … Expand. View 2 excerpts, cites background. In the paper we present a numerical method for solving stiff control problems for delay differential equations based on the method of steps and the differential transformation method DTM.

Spectra and pseudospectra of neutral delay differential equations with application to real-time substructuring. This paper deals with the computation of pseudospectra of neutral delay differential equations NDDEs with fixed finite delays. This method provides information on the sensitivity of eigenvalues of … Expand.

A spectral element approach for the stability analysis of time-periodic delay equations with multiple delays. Nonlinear Sci. Legendre wavelet solution of high order nonlinear ordinary delay differential equations.

The purpose of this paper is to illustrate the use of the Legendre wavelet method in the solution of high-order nonlinear ordinary differential equations with variable and proportional delays. The … Expand. Numerical methods for delay differential equations. The main purpose of the book is to introduce the readers to the numerical integration of the Cauchy problem for delay differential equations DDEs.

Peculiarities and differences that DDEs exhibit … Expand. View 2 excerpts, references background and methods. View 2 excerpts, references methods and background. We first consider the numerical integration of ordinary differential equations ODEs with Runge-Kutta methods that have continuous extensions.

For some methods of this kind we develop robust and … Expand. View 3 excerpts, references background and methods. Software for the numerical solution of systems of functional differential equations with state-dependent delays. Abstract The theoretical basis for the numerical solution of a general class of functional differential equations is reviewed.

A software package for the solution of differential equations with … Expand. View 1 excerpt, references background. Consider again the PDE of Eq. The boundary points, with known values, are excluded. This system of equations is a tridiagonal system, since each equation has three consecutive nonzeros centered around the diagonal.

An approach which requires the solution of simultaneous equations is called an implicit algorithm. A sketch of the C-N stencil is shown in Fig. Note that the coefficient matrix of the C-N system of equations does not change from step to step. Thus, one could compute and save the LU factors of the coefficient matrix, and merely do the forward-backward substitution FBS at each new time step, thus speeding up the calculation.

This speedup would be particularly significant in higher dimensional problems, where the coefficient matrix is no longer tridiagonal. It can be shown that the C-N algorithm is stable for any r, although better accuracy results from a smaller r.

A smaller r corresponds to a smaller time step size for a fixed spacial mesh. C-N also gives better accuracy than the explicit approach for the same r. The only difference between this problem and the one considered earlier in Eq. Assume a mesh labeling as shown in Fig.

Note that c has the dimension of velocity. It can be shown by direct substitution into Eq. The two waves [each equal to half the initial shape f x ] travel in opposite directions from each other at speed c. If f x is nonzero only for a small domain, then, after both waves have passed the region of initial disturbance, the string returns to its rest position. The travelling wave G x is given by Eq.

Thus, as time advances, the center section of the string reaches a state of rest, but not in its original position Fig. From the preceding discussion, it is clear that disturbances travel with speed c.

For an observer at some fixed location, initial displacements occurring elsewhere pass by after a finite time has elapsed, and then the string returns to rest in its original position. Nonzero initial velocity disturbances also travel at speed c, but, once having reached some location, will continue to influence the solution from then on. Figure Domains of Influence and Dependence. Conversely, the domain of dependence of the solution on the initial data consists of all points within a distance ct of the solution point.

If the initial data are discontinuous as, for example, in shocks , the most accurate and the most convenient approach for solving the equations is probably the method of characteristics.

On the other hand, problems without discontinuities can probably be solved most conveniently using finite difference and finite element techniques. Here we consider finite differences. A central finite difference approximation to the PDE, Eq. Thus, this algorithm is an explicit algorithm. The corresponding stencil is shown in Fig. That is, the numerical solution would ignore necessary information. The explicit finite difference algorithm is given in Eq.

If we substitute this last result into Eq. Thus, to implement the explicit finite difference algorithm, we use Eq. At such a boundary, a suitable boundary condition must be imposed to ensure that outgoing waves are not reflected.

Consider a vibrating string which extends to infinity for large x. We truncate the compu- tational domain at some finite x. Let the initial velocity be zero.

Note that the x direction is normal to the boundary. The boundary condition, Eq. This condition is exact in 1-D i. The nonreflecting boundary condition, Eq. For example, at the typical point i, j on the nonreflecting boundary in Fig.

The initial displacement is a triangular-shaped pulse in the middle of the string, similar to Fig. Notice that the triangular wave is absorbed without any reflection from the two boundaries. For steady-state wave motion, the solution u x, t is time-harmonic, i. The nonreflecting boundary condition can be interpreted physically as a damper dash- pot. Consider, for example, a bar undergoing longitudinal vibration and terminated on the right end with the nonreflecting boundary condition, Eq.

Thus, from Eq. Since, from Eq. A mechanical device which applies a force propor- tional to velocity is a dashpot. The minus sign in this equation means that the force opposes the direction of motion, as required to be physically realizable.

Thus, the application of this dashpot to the end of a finite length bar simulates exactly a bar of infinite length Fig. This problem corresponds physically to two-dimensional steady-state heat conduction over a rectangular plate for which temperature is specified on the boundary.

We attempt an approximate solution by introducing a uniform rectangular grid over the domain, and let the point i, j denote the point having the ith value of x and the jth value of y Fig.

Then, using central finite difference approximations to the second derivatives Fig. For example, consider the mesh shown in Fig. Although there are 20 mesh points, 14 are on the boundary, where the temperature is known.

Thus, the resulting numerical problem has only six degrees of freedom unknown variables. The application of Eq.

This linear system of six equations in six unknowns can be solved with standard equation solvers. Because the central difference operator is a 5-point operator, systems of equations of this type would have at most five nonzero terms in each equation, regardless of how large the mesh is.

Thus, for large meshes, the system of equations is sparsely populated, so that sparse matrix solution techniques would be applicable.

Since the numbers assigned to the mesh points in Fig. Some equation solvers based on Gaussian elimination oper- ate more efficiently on systems of equations for which the nonzeros in the coefficient matrix are clustered near the main diagonal. Such a matrix system is called banded. Systems of this type can also be solved using an iterative procedure known as relaxation, which uses the following general algorithm: 1.

Initialize the boundary points to their prescribed values, and initialize the interior points to zero or some other convenient value e. Loop systematically through the interior mesh points, setting each interior point to the average of its four neighbors. Continue this process until the solution converges to the desired accuracy.

For example, consider again the problem of the last section but with a Neumann, rather than Dirichlet, boundary condition on the right side Fig. The procedure is used in a variety of applications, including structural mechanics and dy- namics, acoustics, heat transfer, fluid flow, electric and magnetic fields, and electromagnetics.

Although the main theoretical bases for the finite element method are variational principles and the weighted residual method, it is useful to consider discrete systems first to gain some physical insight into some of the procedures. We let u2 and u3 denote the displacements from the equilibrium of the two masses m2 and m3. The stiffnesses of the two springs are k1 and k2. This approach would be very tedious and error-prone for more complex systems involving many springs and masses.

To develop instead a matrix approach, we first isolate one element, as shown in Fig. In the absence of the masses and constraints, this system is shown in Fig.

This matrix corresponds to the unconstrained system. Notice that, from Eqs. Without constraints, K is singular and the solution of the mechanical problem is not unique because of the presence of rigid body modes. K is symmetric. This property is a special case of the Betti reciprocal theorem in mechanics. An off-diagonal term is zero unless the two points are common to the same element. Thus, K is sparse in general and usually banded.

K is singular without enough constraints to eliminate rigid body motion. For spring systems, that have only one DOF at each point, the sum of any matrix column or row is zero. The forces must sum to zero, since the object is in static equilibrium. We summarize the solution procedure for spring systems: 1. Generate the element stiffness matrices. Assemble the system K and F.

Apply constraints. Compute reactions, spring forces, and stresses. However, matrix assembly for a truss structure a structure made of pin-jointed rods differs from that for a collection of springs, since the rod elements are not all colinear e.

Thus the elements must be transformed to a common coordinate system before the element matrices can be assembled into the system stiffness matrix. A typical rod element in 2-D is shown in Fig. The four DOF for the element are also shown in Fig.

These four values complete the first column of the matrix. Later, in the discussion of matrix transformations, we will derive a more convenient way to obtain this matrix. By using matrix partitioning, we can treat nonzero constraints and recover the forces of constraint. The grid point forces Fs at the constrained DOF consist of both the reactions of constraint and the applied loads, if any, at those DOF.

For example, consider the beam shown in Fig. When the applied load is distributed to the grid points, the loads at the two end grid points would include both reactions and a portion of the applied load. However, for general purpose software, such a capability is essential. An alternative approach to enforce this constraint is to attach a large spring k0 between DOF i and ground a fixed point and to apply a force Fi to DOF i equal to k0 U.

This spring should be many orders of magnitude larger than other typical values in the stiffness matrix. A large number placed on the matrix diagonal will have no adverse effects on the conditioning of the matrix. The large spring approach can be used for any system of equations for which one wants to enforce a constraint on a particular unknown.

The main advantage of this approach is that computer coding is easier, since matrix sizes do not have to change. Choose the large spring k0 to be, say, 10, times the largest diagonal entry in the unconstrained K matrix. For each DOF i which is to be constrained zero or not , add k0 to the diagonal entry Kii , and add k0 U to the corresponding right-hand side term Fi , where U is the desired constraint value for DOF i.

In two dimensions, we designate the four DOF associated with flexure as shown in Fig. Rotations are expressed in radians. For this element, note that there is no coupling between axial and transverse behavior. Transverse shear, which was ignored, can be added. The three-dimensional counterpart to this matrix would have six DOF at each grid point: three translations and three rotations ux , uy , uz , Rx , Ry , Rz.

In addition, for bending in two different planes, there would have to be two moments of inertia I1 and I2 , in addition to a torsional constant J and possibly a product of inertia I For most problems of interest in engineering, exact stiffness matrices cannot be derived. This figure also shows the domain modeled with several triangular elements. A typical element is shown in Fig.

Note that the number of undetermined coefficients equals the number of DOF 6 in the element. The linear approxima- tion in Eq. At the vertices, the displacements in Eq. Hence it is integrable in the neighbourhood of x0. Integrating the initial value problem 2. Since the equivalence of the initial value problem?? Now using the method of iteration, we build up the approximate solutions of the equation 2.

This is done by taking a crude approximation to a solution and feeding it into the equation to be solved to obtain a better approximation to the exact solution. This procedure is done repeatedly as shown in [RCF02]. The proof of this theorem is broken into five steps. From equation 2. To obtain equation 2. Hence by the above, we have proved that equation 2. Also since f is a continuous function in the closed and bounded region R then it is also uniformly continuous in the same region.

By equation 2. Assuming that equation 2. Hence this establishes equation 2. Therefore we have shown that equation 2. Thus, the question of the initial value problem 2. In the above equation 2. Thus our main aim is to show that the unique solution of the initial value problem depends continuously on the initial value y0 [HS99].

So, in this section, we will discuss the continuous dependence of solution on data as presented by [HS99]. So, the question of the initial value problem 2. Maximum Interval Existence of Solutions Having specified through the choice of the real number h which is the interval on which a unique solution exists, and that the true interval of existence of the unique solution is usually larger, we are motivated by the question: what is the maximal interval on which the solution can be extended [HS99]?.

Therefore, we have clearly shown that the initial value problem of the first-order differential equation has a unique solution which depends continuously on the initial conditions in a given maximum interval in which the solution can be extended. Numerical Methods for solving ordinary differential equations Numerical methods have their roots from the Taylor Series. If the differential equation is presented to a computer using a given procedure and form, the computer program will perform the evaluation for the given values.

With the availability of such facilities, we can come up with many methods which have their roots from the Taylor Series[Kre10]. Thus, by successive differentiation, functions such as f1 x, y , f2 x, y ,. Suppose the step size h has a uniform value for all n , then f xn , yn is denoted by fn , simplifying the equation 3.

In each step, we obtain results and use this result to execute the next step that follows. By doing this, we get a sequence of values y0 , y1 , y2 ,. Step 1: define f x, y Step 2: input initial values x0 and y0 Step 3: input h the step size and n the number of steps. On the other hand, the Improved Euler method uses the trapezoidal rule for its approximate values. Section 3. Improved Euler Method Page 14 Considering equation 3. Since Euler method does not give the best and accurate approximate values, we come up with a better method, the Improved Euler method.

This obtained by replacing the integrand of equation 3. Computer algorithm for improved Euler method The following is the computer algorithm for the improved Euler method as shown in [BD12]. Fourth order Runge—Kutta Method Page 15 Step 2: input initial values x0 and y0 Step 3: input h the step size and n the number of steps. The Runge Kutta method originates from the Euler and improved Euler method and is considered to be a member of the family of Runge Kutta methods.

This method is preferred because of its relative simplicity ,sufficiency,accuracy and efficiency when handling a given number of problems. In this discussion, we assume that the step size h is constant [BD12]. Adams—Bashforth Method Page 17 3. Adams—Bashforth Method Page 18 Description of Adams—Bashforth Method Adams—Bashforth methos is more accurate compared to other one-step formulas can be obtained by using a polynomial of a higher degree with more corresponding data points.

Suppose we choose a polynomial P3 x of degree three. Thus, the fourth order formula has a local truncation error proportional to O h5 [BD12].

Adams—Moulton Method Page 19 3. By using an approximating polynomial of a higher degree and correspondingly more data points, we can obtain more accurate higher order formulas. This combination is done so as to obtain accurate and simple results since the pairs are of the same order. Computer algorithm for Predictor-Corrector method The following is the computer algorithm for the predictor-corrector method as shown in [BD12]. Predictor-Corrector Method Page 21 Step 2: input initial values x0 and y0 Step 3: input h the step size and n the number of steps.

Applications of the Numerical Methods on an illustrative example In this chapter, we perform numerical experiments using the different numerical methods on an illustra- tive example of the initial value problem of a first order ordinary differential equation. Let us consider the following differential equation as shown in [BD12].

Example Page 23 The availability of the exact solution easily helps us to determine the accuracy, efficiency and convergence of any numerical method used in this problem. Therefore, this problem will be used throughout the chapter to illustrate and compare the approximate solution of the different numerical methods to the exact solution. Note that there was no round off error during the analysis of the figures used in the tables below.

They are written exactly as they were given as an output by the computer program. Table 4. Example Page 24 Table 4. The approximate solution and the exact solution were combined to generate results that would enable us compare the accuracy and sufficiency of the two using the Euler numerical method.

The obtained results from this numerical experiment are illustrated in the graph 4. Also, the accuracy of the method is not impressive. Hence the Euler method is less accurate since the approximate solution and the exact solution are not the same for the given initial value problem.

Example Page 25 Variation of y with time. This shows that when compared to the Euler method, the improved Euler method is more accurate and efficient since it yields better and substantial results as it also requires much less effort to compute.

It is also clear that the graphs of the different step sizes also vary in terms of convergence of the approximate solution to the exact solution.