Eulers Method

It will be seen that for the three methods Euler, BI and the trapezium method, the same approximation expression is used for the left-hand side of (4.1) but because of points made in questions (2) and (3) above, the methods are very different. 

Consider Fig. 4.1. The curve is the underlying function f(y) that we are trying to find and we have two pieces of information one point on the curve, here at f = 0, the fat point marked in the figure, and the gradient at any point, for example, at the same point, drawn as a tangent. The procedure is now to find a point yi, at a subsequent chosen time t, for example, St, as shown. The picture represents the mathematical problem of finding the solution to the ode (4.1). 

Expressing this mathematically, / is approximated by the simple two-point formula (3.1), written as 

What we have here, with the Euler method, is the definition of the derivative as pertaining to time t (or nSt) and thus f(y t)) or f(y ) on the right-hand side. For our specific example (4.3), this becomes approximately 

When performing simulations of dynamic processes, the engineer must be concerned with the details of numerical integration. In numerical integration what is required is the projection of information known at a particular value of the independent variable, say time, forward to a new or future value of that independent variable. Consider the diflFerential equation 

If we expand a in a Taylor series about a known condition, we have 

Truncating after the first derivative term gives 

Equation (3.1.5) can therefore be used to compute the value a at a new time, t + At, based upon a given value of the function f x, t) which is evaluated at the old time t. Since we know the initial condition at t = 0, we can start this forward -marching scheme and proceed out to any value of t in steps of At. This numerical algorithm of equation (3.1.5) is known as Euler s method. Unfortunately, the local truncation error for this algorithm is large and is of the order of At. To maintain accuracy with a truncation error of this size, the step size needs to be very small. In practice other methods of greater accuracy have been devised in order to allow for increased step sizes. 

The Euler Method can be generalized to handle a set of differential equations 

Consider Fig. 4.1. The curve is the underlying function f(y) that we are trying to find and we have two pieces of information one point on the curve, here at t = 0, the fat point marked in the figure, and the gradient at any 

It turns out that the order of this approximation is also the global error order of the calculation using the Euler method. An alternative way to proceed is to go from the Taylor expansion for y(t + St), as in (3.3), 

---

The method of Ishida et al [84] includes a minimization in the direction in which the path curves, i.e. along (g/ g -g / gj), where g and g are the gradient at the begiiming and the end of an Euler step. This teclmique, called the stabilized Euler method, perfomis much better than the simple Euler method but may become numerically unstable for very small steps. Several other methods, based on higher-order integrators for differential equations, have been proposed [85, 86]. 

Implicit Methods By using different interpolation formulas involving y, it is possible to cferive imphcit integration methods. Implicit methods result in a nonhnear equation to be solved for y so that iterative methods must be used. The backward Euler method is a first-order method. 

To solve the general problem using the backward Euler method, replace the nonlinear differential equation with the nonhuear algebraic equation for one step. 

The basic idea behind the Euler method is to set the change in w per increment of time as... 

Integration is by successive approximation, using essentially Eulers method. For the first trial. 

Simulation by the improved Euler method has shown that a significant radiative heat transfer must be present before reaction zone migration can be demonstrated. 

The simplest method for integrating eq. (14.23) is the Euler method. A series of steps are taken in the direction opposite to the normalized gradient. 

The Euler method, while extremely inaccurate, is also extremely simple. This method is based on the definition of the derivative... 

Discretization error depends on the step size, i.e., if Ax. —> 0, the algorithm would theoretically be exact. The error for Euler method at step N is 0 N(Ax) and total accumulated error is 0 (Ax), that is, it is a first-order method. 

The modified Euler method needs two initial values y and y, and is given by... 

Another improvement on the basic Euler method is to approximate the slope in the middle of the interval by the average of the slopes at the end points, or... 

The value of y is known at some time x, and its value is sought at some future time x + delx. The obvious approach, which is called the direct Euler method, is to approximate the function by a straight line of slope yp(x,y(x)) and to calculate the new value of y from... 

Fig. 2-2. The recovery of atmospheric carbon dioxide calculated by the direct Euler method. The solid line is the analytical solution, and the lines with markers...

Program DGC02 implements this solution, which is called the reverse Euler method. The only change in the code is in the expression foryp in the subroutine EQUATIONS. 

The results are shown in Figure 2-3, in which the solid line is the exact solution. This numerical approach shows no sign of instability even for a time step of 40 years, nearly five times larger than the residence time of atmospheric carbon dioxide (distime). In fact, the reverse Euler method is nearly always stable, and so I shall use it from now on. 

A high degree of accuracy is not called for in many calculations of the evolution of environmental properties because the mathematical description of the environment by a reasonably small number of equations involves an approximation quite independent of any approximation in the equations solution. Figure 2-3 shows how the accuracy of the reverse Euler method degrades as the time step is increased, but it also shows the stability of the method. Even a time step of 40 years, nearly five times larger than the residence time of 8.64 years, yields a solution that behaves like the true solution. In contrast, Figure 2-2 shows the instability of the direct Euler method a time step as small as 10 years introduces oscillations that are not a property of the true solution. 

More complicated numerical methods, such as the Runge-Kutta method, yield more accurate solutions, and for precisely formulated problems requiring accurate solutions these methods are helpful. Examples of such problems are the evolution of planetary orbits or the propagation of seismic waves. But the more accurate numerical methods are much harder to understand and to implement than is the reverse Euler method. In the following chapters, therefore, I shall show the wide range of interesting environmental simulations that are possible with simple numerical methods. 

The key feature of the systems to be considered in this book is that they have short memories that is, the effects of perturbations diminish with the passage of time. In the example of this chapter, the carbon dioxide pressure returns to a value of 1 within a century or two of the perturbation, regardless of the size of the initial perturbation. In this kind of system, computational errors do not grow as the calculation proceeds instead, the system forgets old errors. That is why the reverse Euler method is useful despite its simplicity and limited accuracy. The many properties of the environment that are reasonably stable and predictable can, in principle, be described by equations with just this kind of stability, and these are the properties that can be simulated using the computational methods described in this book. 

The previous chapter showed how the reverse Euler method can be used to solve numerically an ordinary first-order linear differential equation. Most problems in geochemical dynamics involve systems of coupled equations describing related properties of the environment in a number of different reservoirs. In this chapter I shall show how such coupled systems may be treated. I consider first a steady-state situation that yields a system of coupled linear algebraic equations. Such a system can readily be solved by a method called Gaussian elimination and back substitution. I shall present a subroutine, GAUSS, that implements this method. 

The more interesting problems tend to be neither steady state nor linear, and the reverse Euler method can be applied to coupled systems of ordinary differential equations. As it happens, the application requires solving a system of linear algebraic equations, and so subroutine GAUSS can be put to work at once to solve a linear system that evolves in time. The solution of nonlinear systems will be taken up in the next chapter. 

The problem is to calculate new values of the variables, y, at a new time, x 4- delx, given the known values of y at time jc. In the reverse Euler method described in Section 2.4,... 

Program DGC04 solves the time-dependent ocean simulation by the reverse Euler method using Gaussian elimination and back substitution, nrow = 5 the number of equations and unknowns... 

The application of the reverse Euler method of solution to a system of coupled differential equations yields a system of coupled algebraic equations that can be solved by the method of Gaussian elimination and back substitution. In this chapter I demonstrated the solution of simultaneous algebraic equations by means of this method and showed how the solution of algebraic equations can be used to solve the related differential equations. In the process, I presented subroutine GAUSS, the computational engine of all of the programs discussed in the chapters that follow. 

My first attempt to calculate the time history of a geochemical system (Section 2.3) used the obvious approach (the direct Euler method) of evaluating the time derivatives and stepping forward. But it was not sue-... 

Note that prod must be evaluated at the future time x + delx and not at x, because this is the reverse Euler method of calculation, not the direct one. 

As a precautionary note, remember that tvar = x + delx in the application of the reverse Euler method of solution but that tvar — x when the dependent variable is being evaluated for printing or plotting. 

Ordinary differential equations the Euler method Quite commonly, differential equations appear in the form... 

---