Numerical Solutions of Differential Equations

Numerical Solutions of Differential Equations

1

0.5

–1.5

–1

–0.5

0

.5

1

0

–0.5

–1

–1.5

S. Amen · P. Bilokon · A. Brinley Codd · M. Fofaria · T. Shah

Supervised by Prof. Jeff Cash Summer Term, 2004

Abstract Differential equations are an important part of many areas of mathematics, from fluid dynamics to celestial mechanics. They are used by applied mathematicians, physicists and engineers to help in the designing of everything from bridges to ballistic missiles. Except for a few special cases, differential equations cannot normally be solved analytically. Instead, there are many numerical methods which have been developed to provide solutions. This report will investigate several methods, including Euler’s method, Runge-Kutta and the Trapezium rule. Concentrating on Kepler’s equations, we will investigate the accuracy of the methods and find a way to reduce the error to an acceptable tolerance level.

Contents 0 Introduction 0.1 Differential Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.2 Analytical Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.3 Numerical Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4 4 6 6

1 Classification of Differential Equations 1.1 Terminology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Reduction of Linear Differential Equations to a First-order System . . .

7 7 8

2 Numerical Methods 2.1 Forward Euler’s Method . . . . . 2.2 Backward Euler’s Method . . . . 2.3 Trapezium Rule Method . . . . . 2.4 Errors in Numerical Methods . . 2.4.1 Forward Euler’s Method . 2.4.2 Backward Euler’s Method 2.4.3 Trapezium Rule Method . 2.5 Symmetric Methods . . . . . . . 2.6 Extension to Higher Order ODE

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

11 12 13 15 15 16 16 17 18 19

3 Periodic Solutions 21 3.1 Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2 Application of Numerical Methods . . . . . . . . . . . . . . . . . . . . . 22 3.3 Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4 Developing Higher Order Methods 4.1 Developing Formulae . . . . . . . . . . . . . . . 4.2 Solving A Basic Problem . . . . . . . . . . . . . 4.3 Kepler’s Equations . . . . . . . . . . . . . . . . 4.4 Calculating the Earth’s Orbit Around the Sun 4.5 Lobatto Formulae . . . . . . . . . . . . . . . . .

1

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

27 27 29 31 31 32

CONTENTS

CONTENTS

5 Modified Trapeizum Rule 35 5.1 Modified Trapezium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.2 Comparing to other methods . . . . . . . . . . . . . . . . . . . . . . . . 36 6 Solving Stiff Equations 39 6.1 An Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 7 Deriving and Using a Sixth Order Numerical Method 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Deriving Constants . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Deriving A, B and α . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4 Deriving C, D, E and F . . . . . . . . . . . . . . . . . . . . . . . . 7.5 Deriving C, D, E and F . . . . . . . . . . . . . . . . . . . . . . . . 7.6 Usage of the Sixth Order Method . . . . . . . . . . . . . . . . . . . 7.7 The Circle Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.8 Kepler’s Equations for the Movement of the Earth About the Sun 7.9 Kepler’s Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . A Biographies A.1 Sir Isaac Newton (1643-1727) . . . . . A.2 Leonhard Euler (1707-1767) . . . . . . A.3 Carle David Tolmé Runge (1856-1927) A.4 Martin Wilhelm Kutta (1867-1944) . . A.5 Johannes Kepler (1571-1630) . . . . . B Maple Worksheets B.1 feuler.mws . . . . . . . B.1.1 Initialization . B.1.2 Introduction . B.1.3 Implementation B.1.4 Examples . . . B.2 beuler.tex . . . . . . . B.2.1 Initialization . B.2.2 Introduction . B.2.3 Implementation B.2.4 Examples . . . B.3 trap.mws . . . . . . . B.3.1 Initialization . B.3.2 Implementation B.4 circle1.mws . . . . . . B.4.1 Initialization . B.4.2 Introduction .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

2

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . .

. . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . .

43 43 44 44 45 46 47 47 49 50

. . . . .

51 51 52 54 55 55

. . . . . . . . . . . . . . . .

58 59 59 59 59 64 73 73 73 73 77 85 85 85 94 94 94

CONTENTS

B.5

B.6

B.7

B.8

B.9

B.10

CONTENTS

B.4.3 Comparison of Numerical Methods . . . B.4.4 Conclusion . . . . . . . . . . . . . . . . circle2.mws . . . . . . . . . . . . . . . . . . . . B.5.1 Introduction . . . . . . . . . . . . . . . B.5.2 Implementation . . . . . . . . . . . . . . B.5.3 Application to “Circle” Problem . . . . B.5.4 Error Assessment . . . . . . . . . . . . . B.5.5 Conclusion . . . . . . . . . . . . . . . . 4thorder.mws . . . . . . . . . . . . . . . . . . . B.6.1 Initialization . . . . . . . . . . . . . . . B.6.2 Implementation . . . . . . . . . . . . . . 4thorderKepler.mws . . . . . . . . . . . . . . . B.7.1 Initialization . . . . . . . . . . . . . . . B.7.2 Implementation . . . . . . . . . . . . . . modtrap.mws . . . . . . . . . . . . . . . . . . . B.8.1 Initialization . . . . . . . . . . . . . . . B.8.2 Implementation . . . . . . . . . . . . . . stiff.mws . . . . . . . . . . . . . . . . . . . . . . B.9.1 Initialization . . . . . . . . . . . . . . . B.9.2 Background . . . . . . . . . . . . . . . . B.9.3 Implementation . . . . . . . . . . . . . . B.9.4 How to Solve Analytically Using Maple 6thorder.mws . . . . . . . . . . . . . . . . . . . B.10.1 Initialization . . . . . . . . . . . . . . . B.10.2 Implementation . . . . . . . . . . . . . .

3

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .

96 139 140 140 141 144 153 156 157 157 157 177 177 177 191 191 191 213 213 213 214 224 227 227 227

Chapter 0

Introduction Calculus is the study of differential and integral equations and their solutions. It provides a technique for modelling real-world systems in a mathematical manner, by using differential and integral equations that relate the state of a system to states that neighbour it in space and time. It was the work of Leibniz and Newton in developing calculus during the 17th century that allowed for these techniques to be invented.

0.1

Differential Equations

An ordinary differential equation is an equation involving a function y of x and its derivatives, y , y , ..., y (n) . A linear ordinary differential equation of order n takes the following form: an (x)

dn y dn−1 y dy + a0 (x)y = f (x) + an−1 (x) n−1 + . . . + a1 (x) n dx dx dx

(0-1)

Some of the problems and questions that can be answered by modelling real-world problems as a system of differential equations include: • Weather Prediction. • Predicting comet behaviour. • Directing rockets into space. • Measuring rate of radioactive decay. Many radioactive materials disintegrate at a rate proportional to the amount present. • What will the population of a certain country be at some arbitrary time in the future? 4

CHAPTER 0. INTRODUCTION

There are several models that have been used to study population dynamics. One common model is the exponential model, in which the rate of change of the population is proportional to the existing population. • Newton’s Law of cooling. This law has been determined from experimental observations and states that the surface temperature of an object changes at a rate proportional to its relative temperature. These situations and others can be modelled by expressing the physical laws that apply to the situation as a set of equations relating derivatives and integrals. These equations commonly contain variables corresponding to time, velocity, acceleration, mass, density, momentum, energy and the basic laws that apply to many of these situations include conservation of mass and energy, and the balance of forces. These equations can then be solved for some of the variables by assuming other variables to be known. This allows us to extract information about the physical situation being modelled. In some cases the solution to a single differential equation or system of differential equation is exact, e.g. the solution could be a power series or trigonometric function. These solutions, if they exist, are very useful. However the group of differential equations having exact known solutions is small and different techniques are required to solve systems of equations for which an exact solution cannot be found. If an exact solution cannot be found, the alternative is to find an approximation of the solution. The range of techniques used to estimate numerical solutions to differential equations is the main focus of this project. A simple example of a differential equation is y = −y

(0-2)

With initial values y(0) = 1, y (0) = 0, the solution is clearly y = cos x. We will be using this later on to demonstrate the accuracy of some of the numerical methods of obtaining solutions. A more complex example we will also be using is Kepler’s Equations of planetary motion. These are y (0-3) y = − 3 (y 2 + z 2 ) 2 z z = − (0-4) 3 2 (y + z 2 ) 2 with initial values y(0) = 1 − e, z(0) = 0, y (0) = 0, z (0) = 1+e 1−e The solution to this system is an ellipse and e denotes eccentricity. As well as linear ordinary differential equations, we can have non-linear differential equations where the derivatives may be raised to some power, and partial differential 5

CHAPTER 0. INTRODUCTION

equations where y may be a function of several variables. We shall look at the nonlinear Kepler’s Equations later. These more complicated equations can then be used for providing mathematical models of the behavior of sound or light waves or describing how the distribution of heat in an object changes over time. In this study we will restrict ourselves to ordinary differential equations.

0.2

Analytical Solutions

In general, it is not possible to solve a differential equation analytically. However there are some special cases for relatively simple equations where certain tricks can be used to arrive at a solution. The most obvious is the case of the simplest differential equation dy = f (x) dx where f (x) is an integrable function. The solution is clearly y = f (x)dx + c

(0-5)

(0-6)

and if the sample values are given for x and y, c can be calculated easily. More complicated equations can sometimes be rearranged into a form which can be integrated straightforwardly, or sometimes tricks such as the ‘integrating factor’ can be used.

0.3

Numerical Solutions

The remaining differential equations which cannot be solved analytically must be solved numerically. In general, this is done by using initial values given for y (for first order ODEs) and y at x = 0 (for second order ODEs) and then estimating the value of y at x = h where h is some small value. The starting value for y is called y0 , the next calculated value is y1 and so on, so the aim is to use a formula for yn+1 in terms of yn . The main methods we will be using are Euler’s method (both forward and backward) and the trapezium rule. In order to analyse how accurate these methods are, we will use the first example given above (equation 0-2) whose solution is known to be y = cos x.

6

Chapter 1

Classification of Differential Equations In this chapter we introduce the terminology and notation, consider the basic classification of differential equations and develop a method for reducing linear differential equations to a first-order system. 1.1

Terminology 2

dy d y If all the differential coefficients in a differential equation are ordinary (e.g. dx , dx2 , etc.) then the differential equation is an ordinary one. If partial differential coefficients ∂2u are involved (such as ∂u ∂x , ∂y 2 , etc.) then it is a partial differential equation. The order of a differential equation is the order of the highest differential coefficient that it contains. Differential equations of order greater than one are called higher-order differential equations. The degree of a differential equation is the power to which the highest order differential coefficient is raised when the equation is rationalized, that is, converted to a form with no fractional powers. A differential equation of order n dy d2 y , dx2 , . . .). is linear if it is linear in the dependent variable y and its derivatives ( dx Otherwise the equation is called non-linear. In general, an ordinary differential equation is a function of y (the dependent variable), x (the independent variable) and the derivatives of y: dy d2 y d3 y (1-1) f x, y, , 2 , 3 , . . . . dx dx dx

Differential equations play an extremely important role in science and engineering. Any problem involving a “rate of change” of a certain value y with respect to another value, say x, can be formulated using differential equations. Typically the resolution of such a problem involves solving a single linear differential equation, or a system of 7

CHAPTER 1. CLASSIFICATION OF DIFFERENTIAL EQUATIONS

linear differential equations; that is, a set of two or more differential equations with an equal number of unknown functions. As we are often concerned with rates of change w.r.t. time (e.g. speed, velocity, nuclear disintegrations per second, etc.), we may often refer to x as time, although in general it does not have to represent this specific physical quantity. Therefore the graph of y against x is often termed the time series plot (e.g. displacement-time graphs), whereas the graph of y against y is often termed the phase plane plot (e.g. velocity-displacement graphs).

1.2

Reduction of Linear Differential Equations to a Firstorder System

Every initial-value problem of the form bn (x)

dn y dn−1 y + b (x) + . . . + b1 (x)y (x) + b0 (x)y = g(x), n−1 dxn dxn−1

(1-2)

with initial conditions y(x0 ) = c0 , y (x0 ) = c1 , . . . , y (n−1) (x0 ) = cn−1 where bn (x) = 0 can be reduced to the first-order matrix system y(x) ˙ = A(x)y(x) + f (x), y(x0 ) = c,

(1-3) (1-4)

where A(x), f (x), c, and the initial x0 are known. The method of reduction is as follows. Step 1. Hence

Rewrite Equation 1.2 so that

dn y dxn

appears on the left-hand side by itself.

dn y dn−1 y = a + . . . + a1 (x)y (x) + a0 (x)y = f (x), n−1 dxn dxn−1 with ai (x) = −bi (x)/bn (x) for 0 ≤ i ≤ n − 1 and f (x) = g(x)/bn (x).

Step 2.

(1-5)

Define n new variables, y1 (x), y2 (x), . . . , yn (x), by the equations

y1 (x) = y(x), y2 (x) =

dy(x) d2 y(x) dn−1 y(x) , y3 (x) = , . . . , y (x) = . n dx dx2 dxn−1

(1-6)

For these variables we have y˙1 (x) = y2 (x)

(1-7)

y˙2 (x) = y3 (x) .. .

(1-8)

y˙ n−1 (x) = yn (x).

8

(1-9) (1-10)


Step 3.

Express dy n /dx in terms of the new variables. We have dn y(x) d dn−1 y(x) = . y˙ n (x) = dx dxn−1 dxn

(1-11)

Hence, by Equation 1.2, dn−1 y(x) + . . . + a1 (x)y(x) ˙ + a0 (x)y(x) + f (x) dxn−1 = an−1 (x)yn (x) + . . . + a1 (x)y2 (x) + a0 (x)y1 (x) + f (x).

y˙ n (x) = an−1 (x)

(1-12) (1-13)

Finally, reorder the terms: y˙ n (x) = a0 (x)y1 (x) + a1 (x)y2 (x) + . . . + an−1 (x)yn (x) + f (x).

(1-14)

Step 4. Equations 1.2 and 1.2 are a system of first-order linear differential equations in y1 (x), y2 (x), . . . , yn (x). This system is equivalent to the single matrix equation y(x) ˙ = A(x)y(x) + f (x) where ⎞ ⎛ y1 (x) ⎜ y2 (x) ⎟ ⎟ ⎜ (1-15) y(x) = ⎜ . ⎟ , ⎝ .. ⎠

⎛

yn (x) ⎛ ⎞ 0 ⎜ 0 ⎟ ⎜ ⎟ ⎜ ⎟ f (x) = ⎜ ... ⎟ , ⎜ ⎟ ⎝ 0 ⎠ f (x) 0 0 0 .. .

1 0 0 .. .

0 1 0 .. .

0 0 1 .. .

(1-16)

··· ··· ···

⎜ ⎜ ⎜ ⎜ y(x) = ⎜ ⎜ ⎜ ⎝ 0 0 0 0 ··· a0 (x) a1 (x) a2 (x) a3 (x) · · · Step 5.

Finally, we need to take care of the ⎛ c0 ⎜ c1 ⎜ c=⎜ . ⎝ ..

⎞

⎟ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ 1 ⎠ an−1 (x)

(1-17)

initial conditions. Define ⎞

cn−1 9

0 0 0 .. .

⎟ ⎟ ⎟, ⎠

(1-18)


Then the initial conditions can be given by the matrix (vector) equation y(x0 ) = c, and we are done. It is worth noting that it is trivial to adapt this procedure for reducing systems of higher-order equations, rather than a single equation.

10

Chapter 2

Numerical Methods There exist many numerical methods for solving differential equations. They differ in accuracy, performance and applicability. In this chapter we consider some of the simpler numerical methods, namely Forward Euler’s method, Backward Euler’s method and the Trapezium Rule. Many differential equations cannot be solved algebraically (analytically). This means that the solution cannot be expressed as the sum of a finite number of elementary functions (polynomials, exponentials, trigonometric and hyperbolic functions, etc.). But even when it is possible to find an algebraic solution, we may be faced with a system of thousands of differential equations. Many mathematical models, such as weather-forecasting models or fluid dynamics simulations, result in systems of differential equations of enormous complexity and size, making the algebraic solution, if one exists, useless for most practical purposes. In these cases we resort to numerical approximations to the true solution. Some numerical schemes are very simple, such as the Euler’s method, others are relatively complex. In this chapter we consider some of the simpler numerical methods. In general, a numerical method for solving an initial-value problem is a procedure that produces approximate solutions at particular points using only the operations of addition, subtraction, multiplication, division, and functional evaluations. We begin by deriving some numerical methods for solving first-order ordinary differential equations, i.e. equations of the form dy = f (x, y) dx describing the function y(x) in terms of its derivative.

11

(2-1)

CHAPTER 2. NUMERICAL METHODS

2.1

Forward Euler’s Method

We shall derive a numerical scheme called Forward Euler’s method (or simply Euler’s method), which provides an approximation of the solution to a first order differential equation (Equation 2-1). From a graphical perspective, a solution of first-order differential equation is a function whose graph is consistent with the direction field. In other words, the graph of the solution is tangent to the slope lines at every point. Euler’s method approximates the solution by a sequence of points which follows the slope lines. Given Equation 2-1 and the initial condition y(x0 ) = y0 , Euler’s method approximates the solution at x1 , x2 , . . ., where the consecutive x’s differ by h, so that x1 = x0 + h, x2 = x1 + h, .. . xn = xn−1 + h, .. . The increment h is called time step. It is a parameter of the method, which determines the accuracy of the approximation. Let yi denote the approximated value of the solution at xi . The idea is to approximate the solution at each step by the slope line. Consider Figure 2.1. We start at the initial point (x0 , y0 ). The straight line segment is the tangent line to the solution curve at (x0 , y0 ). Its slope is y (x0 ), which is, by Equation 2-1, (2-2) y (x0 ) = f (x0 , y0 ). Hence Euler’s method gives us y1 = y0 + hf (x0 , y0 ).

(2-3)

y2 = y1 + hf (x1 , y1 ).

(2-4)

yn+1 = yn + hf (xn , yn ),

(2-5)

Similarly, at (x1 , y1 ) we have In general, so we can summarize this numerical method by the equation y(x + h) = y(x) + hy (1) (x).

(2-6)

Figure 2.1 illustrates the derivation for the first two steps. Assuming that the value at point x is correct, the Forward Euler’s method computes the value at point x + h with a local error is proportional to h2 . The Forward Euler’s method always gives overshoots on the original curve. 12


y True solution

(x2, y2)

y2

Slope f(x1, y1) Slope f(x0, y0)

Approximated solution y1

(x1, y1)

y0

(x0, y0)

h

h

x0

x1

x2

x

Figure 2.1: Schematic illustration of Forward Euler’s method

2.2

Backward Euler’s Method

The Backward Euler’s method also truncates the Taylor series after two terms. The difference is that the derivative is evaluated at point x + h instead of at point h. Consider Figure 2.2. We start at the initial point (x0 , y0 ). To obtain the next point, (x1 , y1 ), we take the derivative at x1 (not at x0 !) and extrapolate it at point (x0 , y0 ). By Equation 2-1, y (x1 ) = f (x1 , y1 ). (2-7) Hence Backward Euler’s method gives us y1 = y0 + hf (x1 , y1 ).

(2-8)

y2 = y1 + hf (x2 , y2 ).

(2-9)

yn+1 = yn + hf (xn+1 , yn+1 ),

(2-10)

Similarly, at (x1 , y1 ) we have

In general,

13


y (x2, y2)

y2 True solution

Slope f(x2, y2) Slope f(x1, y1)

Approximated solution

(x1, y1)

y1

y0

(x0, y0)

h

x0

h

x1

x2

x

Figure 2.2: Schematic illustration of Backward Euler’s method so we can summarize this numerical method by the equation y(x + h) = y(x) + hy (1) (x + h).

(2-11)

Assuming that the value at point x is correct, the Backward Euler’s method computes the value at point x + h with a local truncation error that is proportional to h2 . The Backward Euler’s method always gives undershoots on the original curve. Normally we do not know the derivative at point x + h, although we need it to compute the function value at point x + h. In practice this requires a rearrangement of the equation. We call such a numerical scheme an implicit numerical scheme, as opposed to explicit numerical schemes, such as the Forward Euler method, where no such rearrangement is necessary. For most equations implicit schemes are more stable than explicit schemes because of the undershoots.

14


2.3

Trapezium Rule Method

The Trapezium Rule (also known as Trapezoidal Rule) method is a simple average of the Forward Euler and Backward Euler schemes: y(x + h) = y(x) + h

y (1) (x) + y (1) (x + h) . 2

(2-12)

It can be shown that the local truncation error is proportional to h3 .

2.4

Errors in Numerical Methods

It should be remembered that numerical methods are simply approximations of the true solution, and as such they inevitably introduce errors. We can classify the errors into three types: 1. The local error describes the error that occurs after a single time step of the method. 2. The global error describes the cumulative effect of the local errors after many time steps. 3. The roundoff error arises from the finite precision with which computers store and process numbers. We shall now provide a rigorous definition of local and global errors. Let the set {(xn , yn ) : n = 0, . . . , N } be an approximation to the solution y(x) of the initial value problem. The global error at time xn is the difference between the true solution at xn and the approximation yn : En = y(xn ) − xn .

(2-13)

The absolute error in the numerical approximation at x = xn is |En |, that is, the absolute value of the global error at xn . The local error, is the difference between the true solution y(xn+1 ) and the approximation after one time step, assuming that (xn , yn ) is a point on the graph of the true solution (i.e. assuming that yn is exact). A numerical method is of order n, where n is a positive integer, if the method is exact for polynomials of degree n or less. In other words, if the true solution of an initial-value problem is a polynomial of degree n or less, then the approximate solution and the true solution will be identical for a method of order n. In general, the higher the order, the more accurate the method. We shall now consider each method in turn, deriving the error term and determining the order of the method. 15


2.4.1

Forward Euler’s Method

Consider the Forward Euler’s method. The local error at the nth step is given by Enloc = y(xn+1 ) − [xn + hf (xn , yn )].

(2-14)

By Taylor’s theorem1 , 1 y(x0 + h) = y(x0 ) + hy (1) (x0 ) + h2 y (2) (ξ) (2-15) 2 1 = y0 + hf (x0 , y0 ) + h2 y (2) (ξ), 2 where ξ is a value between x0 and x0 + h. Hence the local error in Forward Euler’s method is 1 (2-16) E loc = h2 y (2) (ξ) ∝ h2 . 2 Hence the following fact: The local error in Forward Euler’s method decreases as the square of the step size. Thus, if the step size is cut in half, then the local error at each step is reduced by a factor of four. A rigorous derivation of the global error is beyond the scope of this report. We shall restrict ourselves to a simple heuristic argument. Suppose that y(0) = y0 and we want to find y(X). The number of steps required is approximately X/h. If the local error is roughly constant for 0 ≤ x ≤ X, and E(X) is the global error at x = X, E(X) ∝ number of iterations × E loc X × h2 ∝ h ∝ h.

(2-17)

The following fact follows: The global error in Forward Euler’s method decreases linearly with the step size. We say that Euler’s method is a first-order method, because the global error is proportional to h (i.e. first-order with respect to h).

2.4.2

Backward Euler’s Method

To determine the local error for the Backward Euler’s method we expand y(x) in a Taylor series about x1 : 1 y(x1 − h) = y(x1 ) − hy (1) (x1 ) + h2 y (2) (χ). 2

(2-18)

1 We assume that f is at least twice differentiable. It can be shown that this holds if f is continuously differentiable.

16

CHAPTER 2. NUMERICAL METHODS But y(x1 − h) = y(x0 ) = y0 , y(x1 ) = y(x0 + h), and y (1) (x1 ) = f (x1 , y1 ), so we can rearrange Equation 2-18 to obtain 1 y(x0 + h) = y0 + hf (x1 , y1 ) − h2 y (2) (χ). 2

(2-19)

Hence the local error in Backward Euler’s method is 1 E loc = − h2 y (2) (χ) ∝ h2 . 2

(2-20)

This error is O(h2 ), similarly to the error for Forward Euler’s method. Hence the following fact: The local error in Backward Euler’s method decreases as the square of the step size. By following the argument given in Section 2.4, we deduce that the global error is O(h). The global error in Backward Euler’s method decreases linearly with the step size. Once again, a rigorous proof of this fact is beyond the scope of this report.

2.4.3

Trapezium Rule Method

We shall now proceed to deduce the error term for the Trapezium Rule. Geometrically, this numerical scheme corresponds to a rather crude polynomial approximation (a straight line) between successive points xi and xi+1 = xi + h, and hence can only be accurate for sufficiently small h, although, as we shall see, it is “better” than the approximation produced by means of both Forward and Backward Euler’s method. We can express the Trapezium Rule as follows: xi+1 h f (x)dx ≈ (f (xi ) + f (xi+1 )) , (2-21) 2 xi By Taylor expansion, f (xi+1 ) = f (xi + h) = f (xi ) + hf (xi ) +

h2 (2) f (xi ) + . . . , 2!

(2-22)

and hence

xi+1

xi

h h2 (2) f (x)dx ≈ h f (xi ) + f (xi ) + f (xi ) + . . . 2 4

(2-23)

while f (x) may be in xi < x < xi+1 as f (x) = f (xi ) + (x − xi )f (xi ) + 17

(x − xi )2 (2) f (xi ) + . . . 2!

(2-24)


to arrive at the exact form: xi+1 h h2 (2) f (x)dx = h f (xi ) + f (xi ) + f (xi ) + . . . . 2 6 xi Comparison of these two forms shows that the local truncation error is 1 h3 loc 3 1 − f (2) (xi ) + . . . = − f (2) (xi ) + . . . . E =h 6 4 12

(2-25)

(2-26)

We can now make the following conclusion: The local error in Trapezium Rule method decreases as the cube of the step size. We used a slightly different notation in this proof compared to that in Sections 2.4 and 1. In essence, we are trying to find a numerical approximation to the integral of f (x, y) with respect to x. The exact integral is the desired solution of the ordinary differential equation. By following the argument given in Section 2.4, we deduce that the global error is O(h2 ). The global error in Trapezium Rule method decreases as the square of the step size. We have thus demonstrated that by taking the average of the Forward and Backward Euler’s methods we obtain a new method – Trapezium Rule – which is, as far as accuracy is concerned, “better” than the original methods. But accuracy is not the only concern. We shall proceed to demonstrate that there are other criteria for choosing numerical methods, such as stability and symmetricity.

2.5

Symmetric Methods

Consider the Trapezium Rule. By replacing h by −h in 2-12 we obtain y(x − h) = y(x) − h

y (x) + y (x − h) . 2

(2-27)

y (X + h) + y (X) . 2

(2-28)

Now set X = x − h. Then y(X) = y(X + h) − h Hence

y (X) + y (X + h) , (2-29) 2 which is precisely the formula for Trapezium Rule, although this time we are obtaining the solution at point x rather than x + h. Thus the substitution h → −h does not “alter” the numerical scheme. Such methods are called symmetric. In Chapter 2.6 we shall discover the importance of symmetric y(X + h) = y(X) + h

18


methods. In general, for problems with periodic solutions asymmetric (i.e. nonsymmetric) methods fail because they lead to non-periodic solutions. Heuristically, symmetric methods are “invariant” regardless of the “direction” (the sign of h), and produce periodic solutions. Neither Forward, nor Backward Euler scheme is symmetric. Substitution h → −h in 2-6 (Forward Euler’s), for example, gives us the equation y(x − h) = y(x) − hy (x).

(2-30)

On setting X = x − h and rearranging we obtain y(X + h) = y(X) + hy (X + h),

(2-31)

which is, in fact the Backward Euler’s method, which is, in a sense, a “mirror image” of the Forward Euler. Note that by averaging the two “mirror” methods, Forward and Backward Euler, we obtained a symmetric method, the Trapezium Rule.

2.6

Extension to Higher Order ODE

For most practical purposes the numerical analyst is required to solve a system of differential equations rather than a single one. For instance, in weather forecasting, air and water flow models, the mathematical models in use rely on systems of thousands of ODEs. Numerical methods for solving first-order equations are easily extended to a system of first-order initial-value problems. These methods are also applicable to most higherorder initial value problems, in particular those that can be transformed to a system of first-order differential equations by the reduction process described in Section 1.1. The standard form for a system of two equations is y = f (x, y, z),

(2-32)

z = g(x, y, z), with initial conditions y(x0 ) = y0 and z(x0 ) = z0 . Likewise, for three equations we have y = f (x, y, z, w), z = g(x, y, z, w),

(2-33)

w = r(x, y, z, w), with initial conditions y(x0 ) = y0 , z(x0 ) = z0 , w(x0 ) = w0 . If f (x, y, z) = z in 2-32 then the system represents the second-order IV problem y (2) = g(x, y, y ); y(x0 ) = y0 , y (x0 ) = z0 . 19


If f (x, y, z, w) = z in 2-33 then the system represents the third-order IV problem = r(x, y, z, w); y(x0 ) = y0 , y (x0 ) = z0 , y (2) (x0 ) = w0 . It is straightforward to generalize the numerical methods to systems of equations. For example, Forward Euler’s method for 2-32 is as follows.

y (3)

yn+1 = yn + hyn , zn+1 = zn +

(2-34)

hzn .

Similarly it is possible to generalize the other methods to systems of two or more equations, as well as apply them to higher-order initial value problems.

20

Chapter 3

Periodic Solutions Some initial value problems have periodic solutions. For such problems many of the standard numerical methods produce highly inaccurate results, failing to capture the periodic behaviour. We investigate this phenomenon by considering a simple test equation. As a result, we discover that, in order to produce satisfactory results on periodic problems, the numerical method must be symmetric. 3.1

Problem

Some initial value problems have periodic solutions. Consider, for instance, the equation d2 y = −y dx2

(3-1)

with initial conditions1 y(0) = 1,

(3-2)

y (0) = 0.

(3-3)

It is easy to see that y = cos x is a solution, and cos is a periodic function. The time series plot (y against x) is a cosinusoid, whereas the phase plane plot (y against y) is a circle: We would like to test some of our numerical methods on this equation, but before we can apply them, we need to re-write this problem as a first-order system: y = z,

y(0) = 1,

z = −y, 1

y (0) = 0.

(3-4) (3-5)

The initial conditions don’t really matter. We have chosen them to make the solution y = cos x.

21

CHAPTER 3. PERIODIC SOLUTIONS

1

0.5

–10 –8

–6

–4

–2 0

2

4

6

8

10

–0.5

–1

Figure 3.1: xy plot We now wish to integrate the above equation using various methods. We invoke our implementations of the Forward Euler, Backward Euler, and Trapezium Rule given in Chapter 1.2, namely GeneralForwardEuler, GeneralBackwardEuler, and GeneralTrapezium. In Maple, the problem can be specified as follows: > u := [x,y,z]; > f := [1,z,-1*y]; > ic := [0,1,0];

3.2

Application of Numerical Methods

We apply each method in turn with three different time steps, namely h = 10−1 , 10−2 , and 10−3 , and produce the corresponding yz plots. For example, for Forward Euler method with h = 0.1 we use the following system of Maple commands: > fe_sol1 := GeneralForwardEuler(u, f, ic, 0.1, 10): > plot_yz := []: > for i from 1 by 1 to 101 do > plot_yz := [op(plot_yz), [fe_sol1[i,2], fe_sol1[i,3]]]: > end do: > plot(plot_yz); The complete Maple worksheet for this problem can be found on the Project Website: http://www.doc.ic.ac.uk/~pb401/DE/ Figure 3.3 shows the phase plots that we obtained for each of the methods, for all 22


1

0.5

–1

–0.5

0

0.5

1

–0.5

–1

Figure 3.2: yz plot values of h. The Trapezium Rule produced periodic graphs for all values of h, whereas Forward Euler and Backward Euler methods failed to produce a satisfactory solution. As h gets smaller, the method remains periodic for some time, but eventually spirals outwards (Forward Euler) or inwards (Backward Euler). Figure 3.4 shows the time series plot for h = 10−1 , all methods being shown on a single set of axes. The plots for the true solution and the values obtained using the Trapezium Rule are extremely close, as shown by the graph. The graphs for the Forward and Backward Euler methods deviate further and further from the true solution as the distance from the initial value point, (0, 1), increases.

3.3

Error Analysis

Figure 3.5 shows the error plot for h = 10−1 . For the Forward and Backward Euler methods, the error2 oscillates with x, producing increasing maxima as x increases. We shall now consider the error for each method at x = 2π (or the closest point returned by the procedure). The plot in Figure 3.6 demonstrates that the error is proportional to h, as each graph is a straight line.

2

The error at x is defined as |true solution at x − numerical solution at x|.

23


1

1

0.5

0.5

1

0.5

–1.5

–1

–0.5

0.5

1

0 –1

–0.5

0

0.5

1

–1

–0.5

0

0.5

1

–0.5 –0.5

–0.5

–1

–1

–1.5

(a) Forward Euler, h = 10−1

–1

(b) Forward Euler, h = 10−2

0.8

(c) Forward Euler, h = 10−3

1

1

0.5

0.5

0.6 0.4 0.2 –0.8

–0.6

–0.4

–0.2

0.2

0.4

0.6

0.8

1

0 –1

–0.5

0.5

1

–1

–0.5

0

0.5

1

–0.2 –0.4 –0.5

–0.5

–1

–1

–0.6 –0.8

(d) Backward Euler, h = 10−1

–1

–0.5

(e) Backward Euler, h = 10−2

(f) Backward Euler, h = 10−3

1

1

1

0.5

0.5

0.5

0

0.5

1

–1

–0.5

0.5

1

–1

–0.5

0.5

–0.5

–0.5

–0.5

–1

–1

–1

(g) Trapezium Rule, h = 10−1

(h) Trapezium Rule, h = 10−2

1

(i) Trapezium Rule, h = 10−3

Figure 3.3: Phase plane plots for Forward and Backward Euler and Trapezium Rule, for different values of h

24


1 y 0.5 2

4

x

6

8

10

0 –0.5 –1 –1.5

True solution Forward Euler’s Method Backward Euler’s Method Trapezium Rule Method

Figure 3.4: xy plots for all methods on a single set of axes

0.6 0.5 0.4 Error

0.3 0.2 0.1

0

2

4

6

8

10

x Forward Euler’s Method Backward Euler’s Method Trapezium Rule Method

Figure 3.5: Error plots for all methods on a single set of axes

25


0.35 0.3 0.25 h

0.2 0.15 0.1 0.05 0

0.02

0.04

0.06

0.08

0.1

Error Forward Euler’s Method Backward Euler’s Method Trapezium Rule Method

Figure 3.6: The correlation between the error and h

26

Chapter 4

Developing Higher Order Methods We can further develop methods to solve differential equations numerical to reduce the order of error. We can consider several fourth order methods (ie. those that give a local truncation error proportional to h5 and a global error proportional to h4 , quite similar in their approach), based on Runge Kutta formulae 4.1

Developing Formulae

As before our motivation is to solve a second order system, such as d2 y = −y dx2

(4-1)

with initial conditions y(0) = 1, y (0) = 0.

(4-2)

and to get a numerical answer which is periodic. As previously we need to split this into a set of first order differential equations. y = z,

y(0) = 1,

z = −y,

y (0) = 0.

(4-3) (4-4)

Let us consider, an algorithm, which introduces the concept of a “half-step”. yn+1 − yn =

h yn+1 + 4yn+1/2 + yn 6 27

(4-5)

CHAPTER 4. DEVELOPING HIGHER ORDER METHODS

and its derivative

h yn+1 + 4yn+1/2 + yn 6

(4-6)

1 h yn+1/2 = [yn+1 + yn ] − [yn+1 − yn ] 2 8

(4-7)

yn+1 − yn =

We must define the “half-step”

and its derivative

1 h yn+1/2 = [yn+1 + yn ] − [yn+1 − yn ] (4-8) 2 8 This construction is quite similar Simpson’s Rule for evaluating integrals. We can from Equation 4-5 using Equation 4-8. We choose these coefficients eliminate yn+1/2 since a Taylor’s Series exapansion shows it has an error proportional to h4 . We could have done a Taylor’s exapnsion, on Equation 4-9 = A(yn+1 + yn ) + Bhyn + Chyn+1 yn+1/2

(4-9)

Furthermore we can re-express the “half-step” by replacing yn+1 −yn . The algorithm for the system can be neatly summarized into three equations below. Essentially we must solve 4-10 and 4-12 on each iteration, then compute 4-11. For labeling, this is denoted by method A. This method is symmetric (replacing h with −h at the beginning should result in the same algorithm). This is a symmetric fourth order Mono Implicit Runge Kutta (MIRK) method [8].

h2 [4y + 2yn ] 12 n+1/2 h − yn = + yn yn+1 + 4yn+1/2 yn+1 6 yn+1 = yn + hyn +

1 h2 yn+1/2 = [yn+1 + yn ] − + yn yn+1 + 4yn+1/2 2 48 However, we can increase the accuracy of the “half-step”. − yn ] + Ch2 [yn+1 − yn ] yn+1/2 = A[yn+1 + yn ] − Bh[yn+1

(4-10) (4-11)

(4-12)

(4-13)

Taylor’s Expansion can be used to create an initial expression. We then simplify it by − yn . replacing yn+1 − yn and yn+1 1 5h h2 [yn+1 − yn ] + [yn+1 − yn ] yn+1/2 = [yn+1 + yn ] − 2 32 64

(4-14)

1 h2 −2yn+1 + 12yn+1/2 + 14yn yn+1/2 = yn + hyn + 2 192

(4-15)

28


Using this more accurate representation, for yn+1/2 , we can create a more accurate algorithm. For labeling purposes we should call this method B. Thus our aim now is to compute the numerical solution of periodic problems using the fourth order accurate methods A (consisting of Equations 4-10, 4-11 and 4-12) and B (consisting of Equations 4-10, 4-11 and 4-15). The Maple procedures used are called FourthOrderFirst and FourthOrderFirstAcc for A and B respectively.

4.2

Solving A Basic Problem

We are now ready to solve Equation 4-1 using both fourth order methods, using a step-size h = 0.1. Note that we substitute y = z and y = −y into the formulae. Plotting the xy plane we find that the solutions are very close for both A and B. We also plot the yz phase plane and do not encounter the spiralling problems found with backward Euler and forward Euler.

1

0.5

0

2

4

6

8

10

–0.5

–1

Fourth Order A Fourth Order B

Figure 4.1: Plotting xy plane It is instructive to also create an error plot (Figure 4.3) to uncover how precisely the respective solutions are with respect to the real solution, which is y = cos(x). It clearly demonstrates the increased accuracy of using B’s better approximation for yn+1/2 .

29


1

y’

–1

0.5

0

–0.5

0.5 y

1

–0.5

–1


Figure 4.2: Plotting the phase plane

1e–06

8e–07

Err

6e–07

4e–07

2e–07

0

2

4

6

8

10

x Fourth Order A Fourth Order B

Figure 4.3: Error plot, notice how the error increase linearly with x and this is what we would expect for our numerical methods

30


4.3

Kepler’s Equations

We can express Kepler’s Equations in terms of two ODEs. The solution will give us the orbit of a body in two dimensions. y d2 y =− 2 2 dx (y + z 2 )3/2 z d2 z =− 2 2 dx (y + z 2 )3/2

(4-16) (4-17)

with initial conditions (e = 0.6 in our case) y (0) = 0. √ z(0) = 0, z (0) = 1+e 1−e

y(0) = 1 − e,

(4-18) (4-19)

Constructing an appropriate set of first order ODEs, we have,

y = a y a = − 2 (y + z 2 )3/2 z = b z b = − 2 (y + z 2 )3/2

(4-20) (4-21) (4-22) (4-23)

with initial conditions (e = 0.6 in our case) y (0) = a(0) = 0. √ z(0) = 0, z (0) = b(0) = 1+e 1−e

y(0) = 1 − e,

(4-24) (4-25)

We modify the Maple procedures to create FourthOrderFirstKepler and FourthOrderFirstKeplerAcc, to enable them to handle the Kepler Equation, for A and B respectively. Solving the problem, using h = 0.1 we create a yz plot as in Figure 4.4, which will give us the orbit. Both A and B give us the familiar ellipse shape we associate with orbits.

4.4

Calculating the Earth’s Orbit Around the Sun

We can use a similar version of Kepler’s Equations to find out the orbit of the earth around the sun, which has different initial conditions and is. Further details can be found in most good Physics based mechanics books and also http://bartok.ucsc. edu/peter/115/math_intro/node22.html. x can be considered as the independent variable representing time. 31


–0.8

–0.6

–0.4

–0.2

0.2

z 0.4

0.6

0.8

–0.5 y –1

–1.5

Fourth Order A Kepler Problem Fourth Order B Kepler Problem

Figure 4.4: Kepler’s Equations

d2 y GM y =− 2 2 dx (y + z 2 )3/2 GM z d2 z =− 2 2 dx (y + z 2 )3/2

(4-26) (4-27)

with initial conditions y(0) = 1, y (0) = 2π. z(0) = 0,

z (0)

= −π

(4-28) (4-29)

begin where G is the gravitational constant and M is the mass of the sun. It happens that GM = 4π 2 . Before solving using A and B methods, we must create a set of first order ODEs (very similar in form to the previous Kepler’s Equations). We use a step-size h = 0.01 in (Figure 4.5).

4.5

Lobatto Formulae

We again create another, set of formulae, by taking a different approach to calculating − yn in yn+1/2 we subsitute yn+1 . yn+1/2 so rather than substituting yn+1 This eventually results in the “half-step” being defined by, yn+1/2 = yn +

h h2 yn + + 3yn −yn+1 + 4yn+1/2 2 48 32

(4-30)


0.5 –1

–0.5

0.5

z

1

1.5

0 –0.5

y

–1

–1.5

–2

Fourth Order A Kepler Earth Problem Fourth Order B Kepler Earth Problem

Figure 4.5: The earth’s orbit around the sun, method A is not as accurate again

0.5 –1

–0.5

0.5

z

1

–0.5

y

–1

–1.5

–2

Fourth Order A Kepler Earth Problem Fourth Order B Kepler Earth Problem

Figure 4.6: The earth’s orbit around the sun now with h = 0.01 If we solve Equation 4-1, we can create an error plot (Figure 4.7), to compare it to B. We see that B provides a better approximation for the “half-step”.

33


1e–06

8e–07

Err

6e–07

4e–07

2e–07

0

2

4

6

8

10

x Legend Fourth Order B Fourth Order Lobatto

Figure 4.7: Comparing different approximations for the half-step

34

Chapter 5

Modified Trapeizum Rule Here we examine the modified trapezium rule and use it to solve Kepler’s Equation and the circle problem. We also compare it to fourth order method B (in Chapter 4) 5.1

Modified Trapezium

Constructing the Modified Trapezium rule is relatively straightforward. We begin with Equation 5-1 which is the Trapezium Rule, we then take the derivative to give us Equation 5-2. h + yn ] (5-1) yn+1 − yn = [yn+1 2 h [y + yn ] 2 n+1 By eliminating yn+1 from Equation 5-1, we have, h h h yn + (yn+1 + yn ) + yn yn+1 − yn = 2 2 2 yn+1 − yn =

(5-2)

(5-3)

Thus we now have Equation 5-4 for calculating yn+1 and Equation 5-5 for calculating its derivative. These equations form the basis of the Modified Trapezium rule. yn+1 = yn + hyn + = yn + yn+1

h2 [y + yn ] 4 n+1

h [y + yn ] 2 n+1

35

(5-4)

(5-5)

CHAPTER 5. MODIFIED TRAPEIZUM RULE

5.2

Comparing to other methods

Here we compare the Modified Trapezium rule with fourth order method A, from chapter 4. We find that the Modified Trapezium rule produces larger errors. The two problems we examine are Kepler’s Equations and also the circle problem.

1

y’

–1

0.5

–0.5

0.5 y

1

–0.5

–1

Fourth Order A Mod Trap

Figure 5.1: Step size h = 0.1 plot, phase plot of modified trapezium and fourth order A. We find that a periodic solution is generated.

36


1

0.5

0

2

4

6

8

10

–0.5

–1

Fourth Order A Mod Trap

Figure 5.2: Step size h = 0.1 plotting both solutions generated by both methods

0.006 0.005 0.004 Err 0.003 0.002 0.001

0

2

4

6

8

10

x Fourth Order A Modified Trapezium

Figure 5.3: Step size h = 0.1, comparing the errors from both methods

37


–0.8

–0.6

–0.4

–0.2

0.2

0.4

z

0.6

0.8

0

–0.5

–1

y

–1.5

–2

Fourth Order A Kepler Problem Modified Trapezium Kepler Problem

Figure 5.4: Step size h = 0.1, solving Kepler’s Equations using both the Modified Trapezium rule and fourth order A and plotting yz plane

38

Chapter 6

Solving Stiff Equations Certain ODEs are classified as stiff. Essentially this is where solutions can change rapidly with respect to the independent variable, yet the other solution can change a lot more slowly 6.1

An Example

We shall consider the following stiff equation, dy d2 y + λy = 0 + (λ + 1) dx2 dx

(6-1)

with initial conditions y(0) = 1, y (0) = −1.

(6-2)

Introducing an auxillary variable z we can construct a set of first order ODEs from the problem. dy =z (6-3) dx dz = −λy − (λ + 1)z dx

(6-4)

y(0) = 1,z(0) = y (0) = −1.

(6-5)

Solving Equation 6-1 analytically results in y = e−x

39

(6-6)

CHAPTER 6. SOLVING STIFF EQUATIONS However we also have an unwanted solution e−λx , which also satisfies Equation 6-1. This solution changes very rapidly around x = 0. We aim to demonstrate how under certain circumstances Forward Euler is not stable when solving 6-1. In fact, Forward Euler is only stable for step sizes that satisfy, 2 λ

h
with(DEtools): > with(plots): > with(linalg): Warning, the name changecoords has been redefined Warning, the name adjoint has been redefined Warning, the protected names norm and trace have been redefined and unprotected

B.1.2

Introduction

This worksheet provides several implementations of Forward Euler’s Method and demonstrates how they can be used to solve simple IV (initial value) problems. circle.mws uses these results to analyze the “circle” problem and determine which of the methods (Forward Euler’s, Backward Euler’s and Trapezium Rule) produces the most satisfactory results.

B.1.3

Implementation

We provide multiple implementations of Forward Euler’s Method, which differ in their applicability to problems of different orders, the form of the parameters, etc. In each specific case, the user should choose the most appropriate implementation. Here is a brief overview: • SimpleForwardEuler is a simple, but reasonably efficient implementation for first order IV problems. • ForwardEuler is similar to the above, but takes the parameters in a more natural form, then converts them for its purposes and proceeds as SimpleForwardEuler.

59

APPENDIX B. MAPLE WORKSHEETS

• SimpleForwardEuler2 is similar to SimpleForwardEuler, but suitable for second order problems. • ForwardEuler2 is similar to ForwardEuler, but suitable for second order problems. • GeneralForwardEuler : a less efficient but a lot more general implementation of the method, which is applicable to problems of any order, subject to their reducibility to the required input form. • MapleForwardEuler wrapper around Maple’s built-in implementation to make it compatible with ForwardEuler and ForwardEuler2. SimpleForwardEuler: a simple implementation for first order IV problems SimpleForwardEuler is a simplistic implementation of Forward Euler’s Method for solving first-order initial value problems. It accepts the function F , the initial values of the variables x0 and y0 , the step size h, and the number of steps N as parameters. Please refer to the Examples subsection that follows to see how this procedure is applied. > SimpleForwardEuler := proc(F, x0, y0, h, N) > local i, L, X; > X := evalf([x0, y0]); > L := X; > for i from 1 to N do > X := X + [h, h * F(op(X)) ]; > L := L, X; > end do; > return matrix(N + 1, 2, [L]); > end proc: ForwardEuler: a more convenient implementation for first order IV problems ForwardEuler uses a more convenient set of arguments. ode is the ordinary differential equation, ic is the initial condition, interval is the interval over which the solution is to be found, and N is, as before, the number of steps. Please refer to the Examples subsection that follows to see how this procedure is applied.

60

APPENDIX B. MAPLE WORKSHEETS > > > > > > > > > > >

ForwardEuler := proc(ode, ic, interval, N) local h, i, t, y, F, L, X; t := lhs(interval); y := op(0, lhs(ic)); h := (op(2, rhs(interval)) - op(1, rhs(interval))) / N; F := unapply(subs(y(t) = _y, solve(ode, diff(y(t), t))), (t, _y)); x0 := evalf(op(lhs(ic))); y0 := evalf(rhs(ic)); return SimpleForwardEuler(F, x0, y0, h, N); end proc:

Warning, ‘x0‘ is implicitly declared local to procedure ‘ForwardEuler‘ Warning, ‘y0‘ is implicitly declared local to procedure ‘ForwardEuler‘

SimpleForwardEuler2: a simple implementation for second order IV problems SimpleForwardEuler2 is a simplistic implementation of Forward Euler’s Method for solving second-order initial value problems. It accepts the functions F and G, the initial values of the variables x0 , y0 and z0 , the step size h, and the number of steps N as parameters. Please refer to the Examples subsection that follows to see how this procedure is applied. > SimpleForwardEuler2 := proc(F, G, x0, y0, z0, h, N) > local i, L, X; > X := evalf([x0, y0, z0]); > L := X; > for i from 1 to N do > X := X + [ h, h * F(op(X)), h * G(op(X)) ]; > L := L, X; > end do; > return matrix(N + 1, 3, [L]); > end proc: ForwardEuler2: a more convenient implementation for second order IV problems ForwardEuler2 uses a more convenient set of arguments. ode is the ordinary differential equation, ic is the initial condition, interval is the interval over which the solution is to be found, and N is, as before, the number of steps. Please refer to the Examples subsection that follows to see how this procedure is applied.

61

APPENDIX B. MAPLE WORKSHEETS > > > > > > > > > > > > >

ForwardEuler2 := proc(ode, ic1, ic2, interval, N) local t, y, h, F, G, x0, y0, z0, X, L, i; t := lhs(interval); y := op(0, lhs(ic1)); h := (op(2, rhs(interval)) - op(1, rhs(interval))) / N; F := (x, _y, _z) -> _z; G := unapply(subs(y(t) = _y, diff(y(t), t) = _z, solve(ode, diff(diff(y(t), t), t))), (t, _y, _z)); x0 := evalf(op(lhs(ic1))); y0 := evalf(rhs(ic1)); z0 := evalf(rhs(ic2)); return SimpleForwardEuler2(F, G, x0, y0, z0, h, N); end proc:

GeneralForwardEuler: a more general implementation that works for higher order IV problems This implementation of Forward Euler’s Method technically should work for an ODE of ANY order, subject to its reducibility to standard form, etc. The idea behind this implementation is to produce a set of first order equations for each of the variables that represent the different order derivatives, then to simultaneously solve them to get the new set of values for each of the derivatives and for y, then continue to get a list of points. This procedure is used to investigate the “circle” problem. Please refer to circle.mws for a comprehensive example.

62

APPENDIX B. MAPLE WORKSHEETS > > > > > > > > > > > > > > > > > > > > > > > > > >

GeneralForwardEuler := proc(vars, f, init, h, interval) local makeeqn, ans, solveEqns, assignList, diffValue, newICSet,i,j,k, ic; ic := init; makeeqn := nops(u); ans := []; ans := [op(ans), ic]; while op(1,ic) < interval do solveEqns := []; assignList := []; for j from 1 by 1 to makeeqn do assignList := [op(assignList), op(j, u) = op(j, ic)]; end do; for i from 1 by 1 to makeeqn do diffValue := eval(op(i, f),assignList); solveEqns := [op(solveEqns),op(i, u) = op(i, ic) + h * (diffValue)]; end do; newICSet := solve(convert(solveEqns, set), convert(u, set)); for k from 1 by 1 to makeeqn do ic := subsop(k = eval(op(k, u), newICSet), ic); end do; ans := [op(ans), ic]; end do; return matrix(ans); end proc:

MapleForwardEuler: a wrapper around Maple’s built-in implementation Maple provides its own implementation of Forward Euler’s Method. The following is merely a wrapper, which takes the same parameters as ForwardEuler and returns the results obtained by Maple’s native ODE solver based on Forward Euler’s Method. Please refer to the Examples subsection that follows to see how this procedure is applied. > MapleForwardEuler := proc(ode, ic, interval, N) > local x0, h, x_list; > x0 := op(1, rhs(interval)); > h := (op(2, rhs(interval)) - op(1, rhs(interval))) / N; > x_list := vector(N + 1, i -> x0 + (i - 1) * h); > dsolve({ode, ic}, y(x), type=numeric, > method=classical, stepsize=h, output=x_list); > end proc: Supplementary procedures The following procedure converts the “x list” returned by Maple’s dsolve to a list of points, which can be plotted using plot.

63

APPENDIX B. MAPLE WORKSHEETS > > > > >

XListToPoints := proc(x_list) local points_array; points_array := x_list[2, 1]; return convert(points_array, ’listlist’); end proc:

MatrixToPoints converts the matrix of results returned by ForwardEuler and ForwardEuler2 to a list of points, which can be plotted using plot. > MatrixToPoints := proc(x_list) > return convert(x_list, ’listlist’); > end proc:

B.1.4

Examples

Example 1: First order IV problem Suppose we are trying to solve the ODE d dx

y(x) = 3 y(x) + 5

with the initial condition given by y(0) = 7. We specify the problem as follows: > F := (x, y) -> 3 * y + 5; > x0 := 0; > y0 := 7; > h := 0.1; > N := 10; F := (x, y) → 3 y + 5 x0 := 0 y0 := 7 h := 0.1 N := 10 Then call SimpleForwardEuler to solve the system for us: > sol1 := SimpleForwardEuler(F, x0, y0, h, N);

64

APPENDIX B. MAPLE WORKSHEETS ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ sol1 := ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

⎤ 0. 7. ⎥ 0.1 9.6 ⎥ ⎥ 0.2 12.98 ⎥ ⎥ 0.3 17.374 ⎥ ⎥ 0.4 23.0862 ⎥ ⎥ 0.5 30.51206 ⎥ ⎥ 0.6 40.165678 ⎥ ⎥ 0.7 52.7153814 ⎥ ⎥ 0.8 69.02999582 ⎥ ⎥ 0.9 90.23899457 ⎦ 1.0 117.8106929

Here is how we could use ForwardEuler to achieve the same result. Notice how the parameters are passed in this case: > ode := diff(y(x), x) = 3 * y(x) + 5; > ic := y(0) = 7; > interval := x=0..1; > N := 10; ode :=

d dx

y(x) = 3 y(x) + 5

ic := y(0) = 7 interval := x = 0..1 N := 10 Now let’s call ForwardEuler : > sol2 := ForwardEuler(ode, ic, interval, N); ⎡ 0. 7. ⎢ 0.1000000000 9.600000000 ⎢ ⎢ 0.2000000000 12.98000000 ⎢ ⎢ 0.3000000000 17.37400000 ⎢ ⎢ ⎢ 0.4000000000 23.08620000 ⎢ sol2 := ⎢ 0.5000000000 30.51206000 ⎢ ⎢ 0.6000000000 40.16567800 ⎢ ⎢ 0.7000000000 52.71538140 ⎢ ⎢ 0.8000000000 69.02999582 ⎢ ⎣ 0.9000000000 90.23899457 1.000000000 117.8106929

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

The results are identical to those of SimpleForwardEuler . This is hardly surprising, because we have not made any alterations to the procedure’s logic. Let us check whether Maple comes up with the same solution: 65


>

sol_maple := MapleForwardEuler(ode, ic, interval, N); ⎡ ⎤ [x, y(x)] ⎤ ⎡ ⎢ ⎥ 0 7. ⎢ ⎥ ⎢ ⎢ 1 ⎥ ⎥ ⎢ ⎢ 9.59999999999999964 ⎥ ⎥ ⎢ ⎢ 10 ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ 1 ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ 5 12.9800000000000004 ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ 3 ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ 10 17.3739999999999988 ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ 2 ⎥ ⎥ ⎢ ⎢ ⎥ 23.0861999999999981 ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ 5 ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ 1 ⎢ ⎥ ⎥ ⎢ sol maple := ⎢ ⎢ 30.5120599999999982 ⎥ ⎥ ⎢ ⎢ 2 ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ 3 ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ 40.1656779999999998 ⎢ ⎢ 5 ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ 7 ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ 10 52.7153813999999984 ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ 4 ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ 69.0299958200000248 ⎥ ⎥ ⎢ ⎢ ⎢ ⎢ 5 ⎥ ⎥ ⎢ ⎢ ⎥ ⎥ ⎢ ⎢ 9 ⎥ ⎥ ⎢ ⎢ 90.2389945660000308 ⎥ ⎥ ⎣ ⎣ 10 ⎦ ⎦ 1

117.810692935800034

By comparing our values to those returned by Maple, we confirm the correctness of our implementations. We can use Maple to obtain an analytical solution for our ODE. This is achieved by using dsolve with different parameters: > sol_true := dsolve({ode, ic}); 5 26 (3 x) e sol true := y(x) = − + 3 3 Many differential equations do not have analytical solutions, or they are too difficult to calculate. In this case we are lucky, as our equation is very simple. To plot a solution that was obtained from a symbolic computation, we use plot. However, we must extract the function first and specify whichever parameters it takes. In our case, x is the only parameter. > sol_true_f := unapply(op(2,sol_true),x); 5 26 (3 x) e sol true f := x → − + 3 3 > plot(sol_true_f, 0..1); 66


160 140 120 100 80 60 40 20 0

0.2

0.4

0.6

0.8

1

How different is this solution from those obtained by Euler’s Forward Method? To answer this question, we may wish to plot the two solutions on a single graph. odeplot provides a way of plotting the numerical solutions returned by dsolve, for example: > odeplot(sol_maple, [x, y(x)], 0..1);

67


120 100 80 y 60 40 20 0

0.2

0.4

0.6

0.8

1

x

However, odeplot is a lot less flexible and powerful than plot, and it is difficult to combine different kinds of graphs with this method. On the other hand, we need to convert the results of our numerical solution to a set of points before they can be plotted with plot. This is done using XListToPoints procedure: > sol_maple_pts := XListToPoints(sol_maple); 1 1 , 9.59999999999999964], [ , 12.9800000000000004], 10 5 2 3 [ , 17.3739999999999988], [ , 23.0861999999999981], 10 5 3 1 [ , 30.5120599999999982], [ , 40.1656779999999998], 2 5 4 7 [ , 52.7153813999999984], [ , 69.0299958200000248], 10 5 9 [ , 90.2389945660000308], [1, 117.810692935800034]] 10 sol maple pts := [[0, 7.], [

We can now plot the symbolic solution and the numerical solution on the same set of axes: > plot([sol_true_f, sol_maple_pts], 0..1, labels=["x", "y"], > legend=["true", "Euler’s Forward Method"]);

68


160 140 120 100 y 80 60 40 20 0

0.2

0.4

0.6

x

0.8

1

true Euler’s Forward Method

As you can see, the further from the initial point, the more the numerical solution deviates from the true one. Note that the value of y returned by Euler’s Forward Method is always less than the true value. Example 2: Second order Suppose we are trying to solve the ODE d2 dx2

y(x) = 2 y(x)

d y(0) = 5. We specify the problem with the initial conditions given by y(0) = 7 and dx as follows: > ode := diff(diff(y(x), x), x) = 2 * y(x);

ode :=

d2 dx2

y(x) = 2 y(x)

Here is how we could express it in terms of the functions F and G, so that we can apply SimpleForwardEuler2 to obtain a solution:

69

APPENDIX B. MAPLE WORKSHEETS > > > > > > > > > >

t := x; solve(ode, diff(diff(y(t), t), t)); F := (x, _y, _z) -> _z; G := unapply(subs(y(t) = _y, diff(y(t), t) = _z, solve(ode, diff(diff(y(t), t), t))), (t, _y, _z)); x0 := 0; y0 := 7; z0 := 5; h := 0.1; N := 10; t := x 2 y(x) F := (x, y, z ) → z G := (x, y, z ) → 2 y x0 := 0 y0 := 7 z0 := 5 h := 0.1 N := 10

Then call SimpleForwardEuler2 to solve the system for us: > sol1 := SimpleForwardEuler2(F, G, x0, y0, z0, h, N); ⎡ ⎤ 0. 7. 5. ⎢ 0.1 ⎥ 7.5 6.4 ⎢ ⎥ ⎢ 0.2 ⎥ 8.14 7.90 ⎢ ⎥ ⎢ 0.3 ⎥ 8.930 9.528 ⎢ ⎥ ⎢ ⎥ 9.8828 11.3140 ⎥ ⎢ 0.4 ⎢ ⎥ 11.01420 13.29056 ⎥ sol1 := ⎢ 0.5 ⎢ ⎥ 15.493400 ⎥ ⎢ 0.6 12.343256 ⎢ ⎥ ⎢ 0.7 13.8925960 17.9620512 ⎥ ⎢ ⎥ ⎢ 0.8 15.68880112 20.74057040 ⎥ ⎥ ⎢ ⎣ 0.9 17.76285816 23.87833062 ⎦ 1.0 20.15069122 27.43090225 Here is how we could use ForwardEuler2 to achieve the same result. Notice how the parameters are passed in this case: > ode := diff(diff(y(x), x), x) = 2 * y(x); > ic1 := y(0) = 7; > ic2 := D(y)(0) = 5; > interval := x=0..1; > N := 10; ode :=

d2 dx2

y(x) = 2 y(x) 70


ic1 := y(0) = 7 ic2 := D(y)(0) = 5 interval := x = 0..1 N := 10 Now let’s call ForwardEuler2 : > sol2 := ForwardEuler2(ode, ic1, ic2, interval, N); ⎡ 0. 7. 5. ⎢ 0.1000000000 7.500000000 6.400000000 ⎢ ⎢ 0.2000000000 8.140000000 7.900000000 ⎢ ⎢ 0.3000000000 8.930000000 9.528000000 ⎢ ⎢ ⎢ 0.4000000000 9.882800000 11.31400000 ⎢ sol2 := ⎢ 0.5000000000 11.01420000 13.29056000 ⎢ ⎢ 0.6000000000 12.34325600 15.49340000 ⎢ ⎢ 0.7000000000 13.89259600 17.96205120 ⎢ ⎢ 0.8000000000 15.68880112 20.74057040 ⎢ ⎣ 0.9000000000 17.76285816 23.87833062 1.000000000 20.15069122 27.43090225

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

The results are identical to those of SimpleForwardEuler2. We can use Maple to obtain an analytical solution for our ODE. This is achieved by using dsolve with different parameters: > ic1 := y(0) = 7; > ic2 := D(y)(0) = 5; > sol_true := dsolve({ode, ic1, ic2}); ic1 := y(0) = 7 ic2 := D(y)(0) = 5 √ √ 7 5 2 (√2 x) 7 5 2 (−√2 x) )e )e +( − sol true := y(x) = ( + 2 4 2 4 How different is this solution from those obtained by Euler’s Forward Method? To answer this question, we may wish to plot the two solutions on a single graph, just as we did in the previous example. > sol_true_f := unapply(op(2,sol_true),x); > sol1_pts := MatrixToPoints(delcols(sol1, 3..3)); > plot([sol_true_f, sol1_pts], 0..1, labels=["x", "y"], legend=["true", > "Euler’s Forward Method"]); √ √ 7 5 2 (√2 x) 7 5 2 (−√2 x) sol true f := x → ( + )e )e +( − 2 4 2 4 71


sol1 pts := [[0., 7.], [0.1, 7.5], [0.2, 8.14], [0.3, 8.930], [0.4, 9.8828], [0.5, 11.01420], [0.6, 12.343256], [0.7, 13.8925960], [0.8, 15.68880112], [0.9, 17.76285816], [1.0, 20.15069122]]

22 20 18 16 y 14 12 10 8 0

0.2

0.4

x

0.6

0.8

1


Just as before, the further from the initial point, the more the numerical solution deviates from the true one. Also, the value of y is always greater for the true solution.

72


B.2

beuler.tex

Backward Euler’s Method S. Amen, P.A. Bilokon, A.J. Brinley Codd, M. Fofaria, T. Shah B.2.1

Initialization

The following commands initialize Maple. The “warnings” that appear have no adverse effects on the execution, so can be safely ignored. > restart; > with(DEtools): > with(plots): > with(linalg): Warning, the name changecoords has been redefined Warning, the name adjoint has been redefined Warning, the protected names norm and trace have been redefined and unprotected

B.2.2

Introduction

This worksheet provides several implementations of Backward Euler’s Method and demonstrates how they can be used to solve simple IV (initial value) problems. circle.mws uses these results to analyze the “circle” problem and determine which of the methods (Forward Euler’s, Backward Euler’s and Trapezium Rule) produces the most satisfactory results.

B.2.3

Implementation

We provide multiple implementations of Backward Euler’s Method, which differ in their applicability to problems of different orders, the form of the parameters, etc. In each specific case, the user should choose the most appropriate implementation. Here is a brief overview: • SimpleBackwardEuler : a simple, but reasonably efficient implementation for first order IV problems. • BackwardEuler is similar to the above, but takes the parameters in a more natural form, then converts them for its purposes and proceeds as SimpleEuler.

73


• SimpleBackwardEuler2 is similar to SimpleBackwardEuler but suitable for second order problems. • BackwardEuler2 is similar to Euler, but suitable for second order problems. • GeneralBackwardEuler : a less efficient but a lot more general implementation of the method, which is applicable to problems of any order, subject to their reducibility to the required input form. SimpleBackwardEuler: a simple implementation for first order IV problems SimpleBackwardEuler is a simplistic implementation of Backward Euler’s Method for solving first-order initial value problems. It accepts the function F , the initial values of the variables x0 and y0 , the step size h, and the number of steps N as parameters. Please refer to the Examples subsection that follows to see how this procedure is applied. > SimpleBackwardEuler := proc(F, x0, y0, h, N) > local i, L, X, xn, yn, yn1, yn1t; > # create starting point > X := [x0, y0]; > L := X; > for i from 1 to N do > # get xn and yn from last iteration > xn := op(1, X); > yn := op(2, X); > # solve for yn1 (because it is implicit) > yn1 := solve( yn1t = yn + h*F(xn + h, yn1t), yn1t); > X := [xn + h, yn1]; > L := L, X; > end do; > return matrix(N+1, 2, [L]); > end proc: BackwardEuler: a more convenient implementation for first order IV problems BackwardEuler uses a more convenient set of arguments. ode is the ordinary differential equation, ic is the initial condition, interval is the interval over which the solution is to be found, and N is, as before, the number of steps. Please refer to the Examples subsection that follows to see how this procedure is applied.

74

APPENDIX B. MAPLE WORKSHEETS BackwardEuler := proc(ode, ic, interval, N) local h, i, t, y, F, x0, y0; t := lhs(interval); y := op(0, lhs(ic)); h := (op(2, rhs(interval)) - op(1, rhs(interval))) / N; F := unapply(subs(y(t) = _y, solve(ode, diff(y(t), t))), (t, _y)); > x0 := evalf(op(lhs(ic))); > y0 := evalf(rhs(ic)); > return SimpleBackwardEuler(F, x0, y0, h, N); > end proc: > > > > > >

SimpleBackwardEuler2: a simple implementation for second order IV problems SimpleBackwardEuler2 is a simplistic implementation of Backward Euler’s Method for solving second-order initial value problems. It accepts the functions F and G, the initial values of the variables x0 , y0 and z0 , the step size h, and the number of steps N as parameters. Please refer to the Examples subsection that follows to see how this procedure is applied. > SimpleBackwardEuler2 := proc(F, G, x0, y0, z0, h, N) > local i, L, X, xn, yn, yn1, yn1t, F1, zn, zn1, zn1t, Q; > # create starting point > X := evalf([x0, y0, z0]); > L := X; > for i from 1 to N do > # get xn, yn, zn from last iteration > xn := op(1, X); > yn := op(2, X); > zn := op(3, X); > # solve for yn1 (because they are implicit) > # a tuple of the form [yn1, zn1] is returned > Q := solve({yn1t = yn + h*F(xn + h, yn1t, zn1t), zn1t = zn + > h*G(xn + h, yn1t,zn1t)}, {yn1t, zn1t}); > if (lhs(Q[1]) = yn1t) then > yn1 := rhs(Q[1]); > zn1 := rhs(Q[2]); > else > yn1 := rhs(Q[2]); > zn1 := rhs(Q[1]); > end if; > X := [xn + h, yn1, zn1]; > L := L, X; > end do; > return matrix(N + 1, 3, [L]); > end proc:

75


BackwardEuler2: a more convenient implementation for second order IV problems BackwardEuler2 uses a more convenient set of arguments. ode is the ordinary differential equation, ic is the initial condition, interval is the interval over which the solution is to be found, and N is, as before, the number of steps. Please refer to the Examples subsection that follows to see how this procedure is applied. > BackwardEuler2 := proc(ode, ic1, ic2, interval, N) > local t, y, h, F, G, x0, y0, z0, X, L, i; > t := lhs(interval); > y := op(0, lhs(ic1)); > h := (op(2, rhs(interval)) - op(1, rhs(interval))) / N; > F := (x, _y, _z) -> _z; > G := unapply(subs(y(t) = _y, diff(y(t), t) = _z, solve(ode, > diff(diff(y(t), t), t))), (t, _y, _z)); > x0 := evalf(op(lhs(ic1))); > y0 := evalf(rhs(ic1)); > z0 := evalf(rhs(ic2)); > return SimpleBackwardEuler2(F, G, x0, y0, z0, h, N); > end proc: GeneralBackwardEuler: a more general implementation that works for higher order IV problems This implementation of Backward Euler’s Method technically should work for an ODE of ANY order, subject to its reducibility to standard form, etc. The idea behind this implementation is to produce a set of first order equations for each of the variables that represent the different order derivatives, then to simultaneously solve them to get the new set of values for each of the derivatives and for y, then continue to get a list of points. This procedure is used to investigate the “circle” problem. Please refer to circle.mws for a comprehensive example.

76

APPENDIX B. MAPLE WORKSHEETS > >

GeneralBackwardEuler := proc(vars, f, init, h, interval) local makeeqn, ans, solveEqns, assignList, diffValue, newICSet,

i, > > > > > > > > > > > > > > > > > > > > > > > >

j, k, ic; ic := init; makeeqn := nops(u); ans := []; ans := [op(ans), ic]; while op(1,ic) < interval do solveEqns := []; assignList := []; for j from 1 by 1 to makeeqn do assignList := [op(assignList), op(j, u) = op(j, ic)]; end do; for i from 1 by 1 to makeeqn do diffValue := eval(op(i, f), assignList); solveEqns := [op(solveEqns), op(i, u) = op(i, ic) + h * (op(i, f))]; end do; newICSet := solve(convert(solveEqns, set),convert(u, set)); for k from 1 by 1 to makeeqn do ic := subsop(k = eval(op(k, u), newICSet), ic); end do; ans := [op(ans), ic]; end do; return matrix(ans); end proc:

Supplementary procedures The following procedure converts the “x list” returned by Maple’s dsolve to a list of points, which can be plotted using plot. > XListToPoints := proc(x_list) > local points_array; > points_array := x_list[2, 1]; > return convert(points_array, ’listlist’); > end proc: MatrixToPoints converts the matrix of results returned by BackwardEuler and BackwardEuler2 to a list of points, which can be plotted using plot. > MatrixToPoints := proc(x_list) > return convert(x_list, ’listlist’); > end proc:

B.2.4

Examples

Example 1: First order IV problem Suppose we are trying to solve the ODE 77

APPENDIX B. MAPLE WORKSHEETS d dx

y(x) = 3 y(x) + 5

with the initial condition given by y(0) = 7. We specify the problem as follows: > F := (x, y) -> 3 * y + 5; > x0 := 0; > y0 := 7; > h := 0.1; > N := 10; F := (x, y) → 3 y + 5 x0 := 0 y0 := 7 h := 0.1 N := 10 Then call SimpleBackwardEuler to solve the system for us: > sol1 := SimpleBackwardEuler(F, x0, y0, h, N); ⎡ ⎤ 0 7 ⎢ 0.1 10.71428571 ⎥ ⎢ ⎥ ⎢ 0.2 16.02040816 ⎥ ⎢ ⎥ ⎢ 0.3 23.60058309 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ 0.4 34.42940441 ⎥ ⎢ ⎥ sol1 := ⎢ 0.5 49.89914916 ⎥ ⎢ ⎥ ⎢ 0.6 71.99878451 ⎥ ⎢ ⎥ ⎢ 0.7 103.5696922 ⎥ ⎢ ⎥ ⎢ 0.8 148.6709889 ⎥ ⎢ ⎥ ⎣ 0.9 213.1014127 ⎦ 1.0 305.1448753 Here is how we could use BackwardEuler to achieve the same result. Notice how the parameters are passed in this case: > ode := diff(y(x), x) = 3 * y(x) + 5; > ic := y(0) = 7; > interval := x=0..1; > N := 10; ode :=

d dx

y(x) = 3 y(x) + 5

ic := y(0) = 7 interval := x = 0..1 N := 10 Now let’s call BackwardEuler : 78


>

sol2 := BackwardEuler(ode, ic, interval, N); ⎡ 0. 7. ⎢ 0.1000000000 10.71428571 ⎢ ⎢ 0.2000000000 16.02040816 ⎢ ⎢ 0.3000000000 23.60058309 ⎢ ⎢ ⎢ 0.4000000000 34.42940441 ⎢ sol2 := ⎢ 0.5000000000 49.89914916 ⎢ ⎢ 0.6000000000 71.99878451 ⎢ ⎢ 0.7000000000 103.5696922 ⎢ ⎢ 0.8000000000 148.6709889 ⎢ ⎣ 0.9000000000 213.1014127 1.000000000 305.1448753

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

The results are identical to those of SimpleBackwardEuler. This is hardly surprising, because we have not made any alterations to the procedure’s logic. We can use Maple to obtain an analytical solution for our ODE. This is achieved by using dsolve with different parameters: > sol_true := dsolve({ode, ic}); 5 26 (3 x) e sol true := y(x) = − + 3 3 Many differential equations do not have analytical solutions, or they are too difficult to calculate. In this case we are lucky, as our equation is very simple. To plot a solution that was obtained from a symbolic computation, we use plot. However, we must extract the function first and specify whichever parameters it takes. In our case, x is the only parameter. > sol_true_f := unapply(op(2,sol_true),x); 5 26 (3 x) sol true f := x → − + e 3 3 > plot(sol_true_f, 0..1);

79


160 140 120 100 80 60 40 20 0

0.2

0.4

0.6

0.8

1

How different is this solution from those obtained by Euler’s Forward Method? To answer this question, we may wish to plot the two solutions on a single graph, just as we did in the previous example. > sol_true_f := unapply(op(2,sol_true),x); > sol1_pts := MatrixToPoints(sol1); > plot([sol_true_f, sol1_pts], 0..1, labels=["x", "y"], legend=["true", > "Euler’s Forward Method"]); 5 26 (3 x) sol true f := x → − + e 3 3 sol1 pts := [[0, 7], [0.1, 10.71428571], [0.2, 16.02040816], [0.3, 23.60058309], [0.4, 34.42940441], [0.5, 49.89914916], [0.6, 71.99878451], [0.7, 103.5696922], [0.8, 148.6709889], [0.9, 213.1014127], [1.0, 305.1448753]]

80


300 250 200 y 150 100 50 0

0.2

0.4

0.6

x

0.8

1


As you can see, the further from the initial point, the more the numerical solution deviates from the true one. Note that the value of y returned by Euler’s Backward Method is always greater than the true value. Example 2: Second order Suppose we are trying to solve the ODE d2 dx2

y(x) = 2 y(x)

d y(0) = 5. We specify the problem with the initial conditions given by y(0) = 7 and dx as follows: > ode := diff(diff(y(x), x), x) = 2 * y(x);

ode :=

d2 dx2

y(x) = 2 y(x)

Here is how we could express it in terms of the functions F and G, so that we can apply SimpleBackwardEuler2 to obtain a solution:

81

APPENDIX B. MAPLE WORKSHEETS > > > > > > > > > >

t := x; solve(ode, diff(diff(y(t), t), t)); F := (x, _y, _z) -> _z; G := unapply(subs(y(t) = _y, diff(y(t), t) = _z, solve(ode, diff(diff(y(t), t), t))), (t, _y, _z)); x0 := 0; y0 := 7; z0 := 5; h := 0.1; N := 10; t := x 2 y(x) F := (x, y, z ) → z G := (x, y, z ) → 2 y x0 := 0 y0 := 7 z0 := 5 h := 0.1 N := 10

Then call SimpleBackwardEuler2 to solve the system for us: > sol1 := SimpleBackwardEuler2(F, G, x0, y0, z0, ⎡ 0. 7. 5. ⎢ 0.1 7.653061224 6.530612245 ⎢ ⎢ 0.2 8.475635152 8.225739275 ⎢ ⎢ 0.3 9.487968448 10.12333296 ⎢ ⎢ ⎢ 0.4 10.71459362 12.26625168 ⎢ sol1 := ⎢ 0.5 12.18491713 14.70323511 ⎢ ⎢ 0.6 13.93391902 17.49001891 ⎢ ⎢ 0.7 16.00298052 20.69061501 ⎢ ⎢ 0.8 18.44085921 24.37878685 ⎢ ⎣ 0.9 21.30483459 28.63975377 1.0 24.66205099 33.57216397

h, N); ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

Here is how we could use BackwardEuler2 to achieve the same result. Notice how the parameters are passed in this case: > ode := diff(diff(y(x), x), x) = 2 * y(x); > ic1 := y(0) = 7; > ic2 := D(y)(0) = 5; > interval := x=0..1; > N := 10; ode :=

d2 dx2

y(x) = 2 y(x) 82


ic1 := y(0) = 7 ic2 := D(y)(0) = 5 interval := x = 0..1 N := 10 Now let’s call BackwardEuler2 : > sol2 := BackwardEuler2(ode, ⎡ 0. ⎢ 0.1000000000 ⎢ ⎢ 0.2000000000 ⎢ ⎢ 0.3000000000 ⎢ ⎢ ⎢ 0.4000000000 ⎢ sol2 := ⎢ 0.5000000000 ⎢ ⎢ 0.6000000000 ⎢ ⎢ 0.7000000000 ⎢ ⎢ 0.8000000000 ⎢ ⎣ 0.9000000000 1.000000000

ic1, ic2, interval, N); ⎤ 7. 5. 7.653061224 6.530612245 ⎥ ⎥ 8.475635152 8.225739275 ⎥ ⎥ 9.487968448 10.12333296 ⎥ ⎥ ⎥ 10.71459362 12.26625168 ⎥ ⎥ 12.18491713 14.70323511 ⎥ ⎥ 13.93391902 17.49001891 ⎥ ⎥ 16.00298052 20.69061501 ⎥ ⎥ 18.44085921 24.37878685 ⎥ ⎥ 21.30483459 28.63975377 ⎦ 24.66205099 33.57216397

The results are identical to those of SimpleBackwardEuler2. We can use Maple to obtain an analytical solution for our ODE. This is achieved by using dsolve with different parameters: > ic1 := y(0) = 7; > ic2 := D(y)(0) = 5; > sol_true := dsolve({ode, ic1, ic2}); ic1 := y(0) = 7 ic2 := D(y)(0) = 5 √ √ 7 5 2 (√2 x) 7 5 2 (−√2 x) )e )e +( − sol true := y(x) = ( + 2 4 2 4 How different is this solution from those obtained by Euler’s Backward Method? To answer this question, we may wish to plot the two solutions on a single graph, just as we did in the previous example. > sol_true_f := unapply(op(2,sol_true),x); > sol1_pts := MatrixToPoints(delcols(sol1, 3..3)); > plot([sol_true_f, sol1_pts], 0..1, labels=["x", "y"], legend=["true", > "Euler’s Backward Method"]); √ √ 7 5 2 (√2 x) 7 5 2 (−√2 x) sol true f := x → ( + )e )e +( − 2 4 2 4 83


sol1 pts := [[0., 7.], [0.1, 7.653061224], [0.2, 8.475635152], [0.3, 9.487968448], [0.4, 10.71459362], [0.5, 12.18491713], [0.6, 13.93391902], [0.7, 16.00298052], [0.8, 18.44085921], [0.9, 21.30483459], [1.0, 24.66205099]]

24 22 20 18 y 16 14 12 10 8 0

0.2

0.4

x

0.6

0.8

1

true Euler’s Backward Method

Just as before, the further from the initial point, the more the numerical solution deviates from the true one. Also, the value of y is always smaller for the true solution.

84


B.3

trap.mws

Trapezium Rule S. Amen, P.A. Bilokon, A.J. Brinley Codd, M. Fofaria, T. Shah B.3.1

Initialization

The following commands initialize Maple. The “warnings” that appear have no adverse effects on the execution, so can be ignored. > restart; > with(DEtools): > with(plots): > with(linalg): Warning, the name changecoords has been redefined Warning, the name adjoint has been redefined Warning, the protected names norm and trace have been redefined and unprotected

B.3.2

Implementation

Implementation : Trapezium rule This implementation of the Trapezium rule technically should work for an ODE of ANY order. The idea behind the implementation is produce a set of trapezium equations for each of the variables that represent the different order derivatives, then to simultaneously solve them to get the new set of values for each of the derivatives and for y; then continue to get a list of y points.

85


GeneralTrapezium := proc(u, f, init, h, interval) local makeeqn, ans, solveEqns, assignList, diffValue, newICSet,

i, > > > > > > > > > > > > > > > > > > > > > > > >

j, k, ic; ic := init; makeeqn := nops(u); ans := []; ans := [op(ans),ic]; while op(1,ic) < interval do solveEqns := []; assignList := []; for j from 1 by 1 to makeeqn do assignList := [op(assignList), op(j, u) = op(j, ic)]; end do; for i from 1 by 1 to makeeqn do diffValue := eval(op(i, f),assignList); solveEqns := [op(solveEqns), op(i, u) = op(i, ic) + h / 2 * (op(i, f) + diffValue)]; end do; newICSet := solve(convert(solveEqns, set),convert(u, set)); for k from 1 by 1 to makeeqn do ic := subsop(k = eval(op(k, u), newICSet), ic); end do; ans := [op(ans), ic]; end do; return matrix(ans); end proc:

Now to try the trapezium rule on the following ODE : y’ = z, z’ = -y, y(0) = 1, y’(0) = 0 first to put it into a form admissible into our procedure. The procedure then produced estimates for y and y’. This was performed three times by setting h = 0.1, 0.01, and 0.001 respectively. The plot for h = 0.001 was not a complete circle because the number of calculations that needed to be performed was too great.

> > > >

u := [x,y,z]; f := [1,z,-1*y]; ic := [0,1,0]; soln := GeneralTrapezium(u,f,ic,0.1,10); u := [x, y, z] 86

APPENDIX B. MAPLE WORKSHEETS f := [1, z, −y] ic := [0, 1, 0] 0 0.1000000000 0.2000000000 0.3000000000 0.4000000000 0.5000000000 0.6000000000 0.7000000000 0.8000000000 0.9000000000 1. 1.100000000 1.200000000 1.300000000 1.400000000 1.500000000 1.600000000 1.700000000 1.800000000 1.900000000 2. 2.100000000 2.200000000 2.300000000 2.400000000 2.500000000 2.600000000 2.700000000 2.800000000 2.900000000 3. 3.100000000 3.200000000 3.300000000 3.400000000 3.500000000 3.600000000 3.700000000

1 0.9950124688 0.9800996262 0.9554102288 0.9211905548 0.8777819475 0.8256174105 0.7652172884 0.6971840762 0.6221964094 0.5410022946 0.4544116482 0.3632882173 0.2685409637 0.1711149972 0.07198214798 -0.02786872762 -0.1274416109 -0.2257432562 -0.3217930985 -0.4146330346 -0.5033369803 -0.5870201082 -0.6648476739 -0.7360433426 -0.7998969331 -0.8557715018 -0.9031096964 -0.9414393155 -0.9703780188 -0.9896371409 -0.9990245708 -0.9984466687 -0.9879091984 -0.9675172722 -0.9374743009 -0.8980799651 -0.8497272257

87

0 -0.09975062344 -0.1985062282 -0.2952817209 -0.3891117600 -0.4790603851 -0.5642303530 -0.6437720879 -0.7168921561 -0.7828611804 -0.8410211156 -0.8907918127 -0.9316768060 -0.9632682651 -0.9852510632 -0.9974059205 -0.9996115916 -0.9918460747 -0.9741868314 -0.9468100137 -0.9099887071 -0.8640902064 -0.8095723520 -0.7469789629 -0.6769344121 -0.6001373983 -0.5173539765 -0.4294099166 -0.3371824660 -0.2415915993 -0.1435908414 -0.04415775582 0.05571580615 0.1550335995 0.2528049230 0.3480545016 0.4398322149 0.5272225745


3.800000000 3.900000000 4. 4.100000000 4.200000000 4.300000000 4.400000000 4.500000000 4.600000000 4.700000000 4.800000000 4.900000000 5. 5.100000000 5.200000000 5.300000000 5.400000000 5.500000000 5.600000000 5.700000000 5.800000000 5.900000000 6. 6.100000000 6.200000000 6.300000000 6.400000000 6.500000000 6.600000000 6.700000000 6.800000000 6.900000000 7. 7.100000000 7.200000000 7.300000000 7.400000000 7.500000000 7.600000000 7.700000000 7.800000000

-0.7928984042 -0.7281603717 -0.6561588941 -0.5776121906 -0.4933037695 -0.4040746125 -0.3108147860 -0.2144545627 -0.1159551418 -0.01629906110 0.08351960374 0.1825051553 0.2796702065 0.3740455300 0.4646897261 0.5506986132 0.6312142473 0.7054334798 0.7726159694 0.8320915666 0.8832669986 0.9256317872 0.9587633410 0.9823311707 0.9961001857 0.9999330392 0.9937914983 0.9777368253 0.9519291665 0.9166259549 0.8721793423 0.8190326864 0.7577161284 0.6888413048 0.6130952462 0.5312335243 0.4440727149 0.3524822525 0.2573757577 0.1597019236 0.06043505280

88

0.6093538560 0.6854067948 0.7546227581 0.8163113123 0.8698571103 0.9147260294 0.9504704993 0.9767339667 0.9932544519 0.9998671621 0.9965061348 0.9832048968 0.9600961287 0.9274103419 0.8854735791 0.8347041621 0.7756085190 0.7087761326 0.6348736601 0.5546382833 0.4688703551 0.3784254158 0.2842056594 0.1871509339 0.08822936611 -0.01157229513 -0.1112585220 -0.2098349382 -0.3063182378 -0.3997459938 -0.4891862586 -0.5737468600 -0.6525843007 -0.7249121723 -0.7900089998 -0.8472254383 -0.8959907502 -0.9358184985 -0.9663113990 -0.9871652831 -0.9981721319


7.900000000 8. 8.100000000 8.200000000 8.300000000 8.400000000 8.500000000 8.600000000 8.700000000 8.800000000 8.900000000 9. 9.100000000 9.200000000 9.300000000 9.400000000 9.500000000 9.600000000 9.700000000 9.800000000 9.900000000 10.

-0.03943466140 -0.1389110124 -0.2370017174 -0.3327283155 -0.4251359279 -0.5133027829 -0.5963494106 -0.6734474158 -0.7438277411 -0.8067883383 -0.8617011715 -0.9080184818 -0.9452782511 -0.9731088109 -0.9912325497 -0.9994686820 -0.9977350523 -0.9860489528 -0.9645269535 -0.9333837378 -0.8929299611 -0.8435691523

89

-0.9992221519 -0.9903048682 -0.9715092317 -0.9430227300 -0.9051295178 -0.8582075823 -0.8027249727 -0.7392351314 -0.6683713735 -0.5908405695 -0.5074160940 -0.4189301113 -0.3262652746 -0.2303459215 -0.1321288535 -0.03259379192 0.06726639479 0.1664555950 0.2639843903 0.3588799249 0.4501956099 0.5370205656

APPENDIX B. MAPLE WORKSHEETS > > > > >

plot_yz := []: for i from 1 by 1 to 101 do plot_yz := [op(plot_yz),[soln[i,2],soln[i,3]]]: end do: plot(plot_yz);

1

0.5

–1

–0.5

0

0.5

–0.5

–1

90

1


> > > > > >

soln1 := GeneralTrapezium(u,f,ic,0.01,10): plot_yz1 := []: for i from 1 by 1 to 1001 do plot_yz1 := [op(plot_yz1),[soln1[i,2],soln1[i,3]]]: end do: plot(plot_yz1);

1

0.5

–1

–0.5

0.5 –0.5

–1

91

1


> > > > > >

soln2 := GeneralTrapezium(u,f,ic,0.001,1): plot_yz2 := []: for i from 1 by 1 to 1001 do plot_yz2 := [op(plot_yz2),[soln2[i,2],soln2[i,3]]]: end do: plot(plot_yz2);

0.6

0.7

0.8

0 –0.2

–0.4

–0.6

–0.8

92

0.9

1


Now compare the analytical solution of this problem to the solutions obtained using the trapezium rule. > ode := diff(diff(y(x),x),x)=-y(x); > ic1 := y(0) = 1; ic2 := D(y)(0) = 0; > trueSoln := dsolve({ode, ic1, ic2}); ode :=

d2 dx2

y(x) = −y(x)

ic1 := y(0) = 1 ic2 := D(y)(0) = 0 > > > > > > > >

trueSoln := y(x) = cos(x) trueSoln := y(x) = cos(x); plot_real := []: plot_estimate := []: for i from 1 by 1 to 101 do plot_real := [op(plot_real),[soln[i,1], cos(soln [i,1])]]: plot_estimate := [op(plot_estimate),[soln[i,1], soln[i,2]]]: end do: plot([plot_real, plot_estimate]);

1

0.5

0

2

4

6

8

10

–0.5

–1

The plots for the real values and the values obtained using the trapezium rule are incredibly close, as shown by the graph, which contains two lines, although this is not obvious. 93


B.4

circle1.mws

Comparison of IVP Methods I: Periodicity - Part 1 S. Amen, P.A. Bilokon, A.J. Brinley Codd, M. Fofaria, T. Shah B.4.1

Initialization


B.4.2

Introduction

Consider the differential equation d2 y(x) = −y(x) dx2 with initial conditions y(0) = 1 and y (0) = 0. The solution is readily available using analytical methods, pencil and paper. In fact, it is y(x) = cos(x), and hence

y (x) = −sin(x).

Let’s formulate this problem in Maple. > ode := diff(diff(y(x), x), x) = -y(x); > ic1 := y(0) = 1; > ic2 := D(y)(0) = 0; ode :=

d2 dx2

y(x) = −y(x)

ic1 := y(0) = 1 ic2 := D(y)(0) = 0 94


Note the use of the differential operator D(f ) for specifying the initial value d dx y(0) = 0. Maple can calculate the true solution algebraically. By default, dsolve resorts to algebraic methods. > desol_true := dsolve({ode, ic1, ic2}); desol true := y(x) = cos(x) We have thus verified our solution. We shall use plot to plot the true solution. However, we must extract the function ( f ) first and specify whichever parameters it takes. In our case, x is the only parameter. > f := unapply(op(2,desol_true),x); f := cos And here is our plot: > plot(f, -10..10);

1

0.5

–10 –8

–6

–4

–2 0

2

4

6

8

10

–0.5

–1

We can also use Maple to produce a plot of > plot([f, D(f), -100..100]);

95

d dx

y(x) against y(x):


1

0.5

–1

–0.5

0

0.5

1

–0.5

–1

As you can see, the graph is a circle. The solution is periodic, the graph repeatedly retraces the same path over the range of x. We shall use this example to evaluate the quality of solutions produced by numerical methods. Can they capture this periodic behaviour? By comparing the plots of d dx y(x) against y(x) for the various numerical methods with the true solution (the circle, as shown above), we can reason about the quality of the models that they provide.

B.4.3

Comparison of Numerical Methods

Maple implements the following methods for solving IVPs numerically: • rkf45 - A Fehlberg fourth-fifth order Runge-Kutta method with degree four interpolant. This is the default non-stiff numeric method. • rosenbrock - An Implicit Rosenbrock third-fourth order Runge-Kutta method with degree three interpolant. This is the default method when stiffness is required. 96


• dverk78 - A seventh-eighth order continuous Runge-Kutta method. This method is capable of working in arbitrary precision based on the setting of Digits, and can be used to obtain high accuracy solutions for ODE systems. • taylorseries - A Taylor Series method. This method can be used for high accuracy solutions. It usually takes more time than other methods for low accuracy results. Therefore, it is better to use this method only when a very high degree of accuracy is desired. • classical - May refer to one of the “classical” numerical methods: • foreuler (default) is the forward Euler method specified by the equation: Y [n + 1] = Y [n] + h ∗ f (t[n], Y [n]) • heunform is the Heun formula (also known as the trapezoidal rule, or the improved Euler method), as specified by the equation: Y [n + 1] = Y [n] + (h/2) ∗ (f (t[n], Y [n]) + f (t[n + 1], Y [n + 1])) • impoly is the improved polygon method (also known as the modified Euler method), as specified by the equation: Y [n + 1] = Y [n] + h ∗ (f (t[n] + h/2, Y [n] + (h/2) ∗ f (t[n], Y [n]))) • rk2 is the second-order classical Runge-Kutta method, as specified by: k1 = f (t[n], Y [n]) k2 = f (t[n] + h, Y [n] + h ∗ k1) Y [n + 1] = Y [n] + (h/2) ∗ (k1 + k2) • rk3 is the third-order classical Runge-Kutta method, as specified by: k1 = f (t[n], Y [n]) k2 = f (t[n] + (h/2), Y [n] + (h/2) ∗ k1) k3 = f (t[n] + h, Y [n] + (2 ∗ h/3) ∗ k2) Y [n + 1] = Y [n] + (h/6) ∗ (k1 + 4 ∗ k2 + k3) 97


• rk4 is the fourth-order classical Runge-Kutta method, as specified by: k1 = f (t[n], Y [n]) k2 = f (t[n] + h/2, Y [n] + (h/2) ∗ k1) k3 = f (t[n] + h/2, Y [n] + (h/2) ∗ k2) k4 = f (t[n] + h, Y [n] + h ∗ k3) Y [n + 1] = Y [n] + (h/6) ∗ (k1 + 2 ∗ k2 + 2 ∗ k3 + k4) This is not to be confused with method=rkf45, which uses a Fehlberg fourth-fifth order Runge-Kutta method. • adambash is the Adams-Bashford method (a “predictor” method), as specified by: Y [n + 1] =Y [n] + (h/24) ∗ (55 ∗ f (t[n], Y [n]) − 59 ∗ f (t[n − 1], Y [n − 1]) + 37 ∗ f (t[n − 2], Y [n − 2]) − 9 ∗ f (t[n − 3], Y [n − 3]))

• abmoulton is the Adams-Bashford-Moulton method (a “predictor-corrector” method), as specified by: Y [n + 1] =Y [n] + (h/24) ∗ (9 ∗ f (t[n + 1], Y [n + 1]) + 19 ∗ f (t[n], Y [n]) − 5 ∗ f (t[n − 1], Y [n − 1]) + f (t[n − 2], Y [n − 2])), where f(tn+1 , Yn+1 ) is found by first applying the Adams-Bashford method (the “predictor”), then using the above Adams-Bashford-Moulton method (the “corrector”). Stiffness is difficult to define because it involves both the IVP and the method applied to solve it. A stiff IVP has a solution that is very stable (some solutions that start near it converge to it very quickly) and is slowly varying (other solutions change on much faster time scales). Special methods allow stiff problems to be solved efficiently, but they can be very inefficient for problems that are not stiff. Let us apply each of these methods in turn to our problem: > desol_rkf45 := dsolve({ode, ic1, ic2}, numeric, method=rkf45): > desol_rosenbrock := dsolve({ode, ic1, ic2}, numeric, > method=rosenbrock): > desol_dverk78 := dsolve({ode, ic1, ic2}, numeric, > method=dverk78): 98


> > > > > > > > > > > > > > > > > > >

desol_lsode := dsolve({ode, ic1, ic2}, numeric, method=lsode): desol_taylorseries := dsolve({ode, ic1, ic2}, numeric, method=taylorseries): desol_classical_foreuler := dsolve({ode, ic1, ic2}, numeric, method=classical[foreuler]): desol_classical_heunform := dsolve({ode, ic1, ic2}, numeric, method=classical[heunform]): desol_classical_impoly := dsolve({ode, ic1, ic2}, numeric, method=classical[impoly]): desol_classical_rk2 := dsolve({ode, ic1, ic2}, numeric, method=classical[rk2]): desol_classical_rk3 := dsolve({ode, ic1, ic2}, numeric, method=classical[rk3]): desol_classical_rk4 := dsolve({ode, ic1, ic2}, numeric, method=classical[rk4]): desol_classical_adambash := dsolve({ode, ic1, ic2}, numeric, method=classical[adambash]): desol_classical_abmoulton := dsolve({ode, ic1, ic2}, numeric, method=classical[abmoulton]):

We shall now plot our graphs for the short range of x (-10..10), medium range (-100..100) and long range (-500..500).

99


Short Range of x > odeplot(desol_rkf45, [y(x), diff(y(x), x)], x=-10..10);

1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

100

1


>

odeplot(desol_rosenbrock, [y(x), diff(y(x), x)], x=-10..10);

1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

101

1


>

odeplot(desol_dverk78, [y(x), diff(y(x), x)], x=-10..10);

1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

102

1


>

odeplot(desol_lsode, [y(x), diff(y(x), x)], x=-10..10);

1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

103

1


>

odeplot(desol_taylorseries, [y(x), diff(y(x), x)], x=-10..10);

1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

104

1


odeplot(desol_classical_foreuler, [y(x), diff(y(x), x)], x=-10..10);

1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

105

1


odeplot(desol_classical_heunform, [y(x), diff(y(x), x)], x=-10..10);

1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

106

1

APPENDIX B. MAPLE WORKSHEETS > odeplot(desol_classical_impoly, [y(x), diff(y(x), x)], x=-10..10);

1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

107

1


>

odeplot(desol_classical_rk2, [y(x), diff(y(x), x)], x=-10..10);

1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

108

1


>


1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

109

1


>


1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

110

1


odeplot(desol_classical_adambash, [y(x), diff(y(x), x)], x=-10..10);

1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

111

1


odeplot(desol_classical_abmoulton, [y(x), diff(y(x), x)], x=-10..10);

1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

112

1


Medium Range of x > odeplot(desol_rkf45, [y(x), diff(y(x), x)], x=-100..100);

1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

113

1


>


1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

114

1


>


1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

115

1


>


1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

116

1


>


1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

117

1



1 y’ 0.5

–1

–0.5

0

0.5

1 y

–0.5 –1

118



1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

119

1


odeplot(desol_classical_impoly, [y(x), diff(y(x), x)], x=-100..100);

1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

120

1


>


1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

121

1


>


1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

122

1


>


1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

123

1



1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

124

1



1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

125

1


Long Range of x > odeplot(desol_rkf45, [y(x), diff(y(x), x)], x=-500..500);

1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

126

1


>


1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

127

1


>


1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

128

1


>


1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

129

1


>


1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

130

1



3 2 y’ 1 –3

–2

–1

1

2 y

–1 –2 –3

131

3



1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

132

1


odeplot(desol_classical_impoly, [y(x), diff(y(x), x)], x=-500..500);

1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

133

1


>


1

y’

–1

–0.5

0.5

0

0.5 y

–0.5

–1

134

1


>


1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

135

1


>


1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

136

1



1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

137

1



1

y’

–1

0.5

–0.5

0.5 y –0.5

–1

138

1


B.4.4

Conclusion

For a small range of x all methods seem to be appropriate. However, as x gets further from the initial value, Euler’s forward method (classical) “spirals out” producing highly inaccurate results. Other methods seem to cope well with this situation. Please note that we have considered a single aspect of the problem, accuracy, over a limited range of x. We did not assess any other parameters, including the performance of the methods, which distinguish them from each other. Further work is certainly required.

139


B.5

circle2.mws

Comparison of IVP Methods I: Periodicity - Part 2 S. Amen, P.A. Bilokon, A.J. Brinley Codd, M. Fofaria, T. Shah B.5.1

Introduction

Consider the differential equation d2 y(x) = −y(x) dx2 with initial conditions y(0) = 1 and y (0) = 0. The solution is readily available using analytical methods, pencil and paper. In fact, it is y(x) = cos(x), and hence

y (x) = −sin(x).

We now wish to integrate the above equation using Forward Euler’s Method, Backward Euler’s Method and the Trapezium Rule. Our goal is to determine the most appropriate method from this set. By “most appropriate” we mean, in this case, the method that produces the smallest error and is thus closest to the true solution.

140


B.5.2

Implementation

The following procedures have been copied from feuler.mws, beuler.mws and trap.mws. They implement the Forward Euler’s Method, Backward Euler’s Method and the Trapezium Rule Method, respectively. We shall use the most general implementations, which can be applied to the “circle” problem. In this case we are indifferent to their performance. Forward Euler’s Method > GeneralForwardEuler := proc(vars, f, init, h, interval) > local makeeqn, ans, solveEqns, assignList, diffValue, > newICSet,i,j,k, ic; > ic := init; > makeeqn := nops(u); > ans := []; > ans := [op(ans), ic]; > while op(1,ic) < interval do > solveEqns := []; > assignList := []; > for j from 1 by 1 to makeeqn do > assignList := [op(assignList), op(j, u) = op(j, ic)]; > end do; > for i from 1 by 1 to makeeqn do > diffValue := eval(op(i, f),assignList); > solveEqns := [op(solveEqns),op(i, u) = op(i, ic) + h * > (diffValue)]; > end do; > newICSet := solve(convert(solveEqns, set), convert(u, set)); > for k from 1 by 1 to makeeqn do > ic := subsop(k = eval(op(k, u), newICSet), ic); > end do; > ans := [op(ans), ic]; > end do; > return matrix(ans); > end proc:

141


Backward Euler’s Method > GeneralBackwardEuler := proc(vars, f, init, h, interval) > local makeeqn, ans, solveEqns, assignList, diffValue, newICSet, i, > j, k, ic; > ic := init; > makeeqn := nops(u); > ans := []; > ans := [op(ans), ic]; > while op(1,ic) < interval do > solveEqns := []; > assignList := []; > for j from 1 by 1 to makeeqn do > assignList := [op(assignList), op(j, u) = op(j, ic)]; > end do; > for i from 1 by 1 to makeeqn do > diffValue := eval(op(i, f), assignList); > solveEqns := [op(solveEqns), op(i, u) = op(i, ic) + h * (op(i, > f))]; > end do; > newICSet := solve(convert(solveEqns, set),convert(u, set)); > for k from 1 by 1 to makeeqn do > ic := subsop(k = eval(op(k, u), newICSet), ic); > end do; > ans := [op(ans), ic]; > end do; > return matrix(ans); > end proc:

142


Trapezium Rule Method > GeneralTrapezium := proc(vars, f, init, h, interval) > local makeeqn, ans, solveEqns, assignList, diffValue, newICSet, i, > j, k, ic; > ic := init; > makeeqn := nops(u); > ans := []; > ans := [op(ans),ic]; > while op(1,ic) < interval do > solveEqns := []; > assignList := []; > for j from 1 by 1 to makeeqn do > assignList := [op(assignList), op(j, u) = op(j, ic)]; > end do; > for i from 1 by 1 to makeeqn do > diffValue := eval(op(i, f), assignList); > solveEqns := [op(solveEqns), op(i, u) = op(i, ic) + h / 2 * > (op(i, f) + diffValue)]; > end do; > newICSet := solve(convert(solveEqns, set),convert(u, set)); > for k from 1 by 1 to makeeqn do > ic := subsop(k = eval(op(k, u), newICSet), ic); > end do; > ans := [op(ans), ic]; > end do; > return matrix(ans); > end proc:

143


B.5.3

Application to “Circle” Problem

Now let’s try the three methods on the “circle” problem. First to put it into a form admissible into our procedure: > u := [x,y,z]; > f := [1,z,-1*y]; > ic := [0,1,0]; u := [x, y, z] f := [1, z, −y] ic := [0, 1, 0] The procedure then produces estimates for y and y’. This was performed three times by setting h = 0.1, 0.01, and 0.001 respectively. The plot for h = 0.001 was not a complete circle because the number of calculations that needed to be performed was too great. FE, h = 0.1 > fe_sol1 := GeneralForwardEuler(u, f, ic, 0.1, 10): > plot_yz := []: > for i from 1 by 1 to 101 do > plot_yz := [op(plot_yz), [fe_sol1[i,2], fe_sol1[i,3]]]: > end do: > plot(plot_yz);

1 0.5 –1.5

–1

–0.5

0.5 0 –0.5 –1 –1.5

144

1


FE, h = 0.01 > fe_sol2 := GeneralForwardEuler(u, f, ic, 0.01, 10): > plot_yz := []: > for i from 1 by 1 to 1001 do > plot_yz := [op(plot_yz), [fe_sol2[i,2], fe_sol2[i,3]]]: > end do: > plot(plot_yz);

1

0.5

–1

–0.5

0

0.5

–0.5

–1

145

1


FE, h = 0.001 > fe_sol3 := GeneralForwardEuler(u, f, ic, 0.001, 10): > plot_yz := []: > for i from 1 by 1 to 1001 do > plot_yz := [op(plot_yz), [fe_sol3[i,2], fe_sol3[i,3]]]: > end do: > plot(plot_yz);

0.6

0.7

0.8

0 –0.2

–0.4

–0.6

–0.8

146

0.9

1


BE, h = 0.1 > be_sol1 := GeneralBackwardEuler(u, f, ic, 0.1, 10): > plot_yz := []: > for i from 1 by 1 to 101 do > plot_yz := [op(plot_yz), [be_sol1[i,2], be_sol1[i,3]]]: > end do: > plot(plot_yz);

0.8 0.6 0.4 0.2 –0.8 –0.6 –0.4 –0.2 0 –0.2

0.2 0.4 0.6 0.8

–0.4 –0.6 –0.8

147

1


FE, h = 0.01 > be_sol2 := GeneralBackwardEuler(u, f, ic, 0.01, 10): > plot_yz := []: > for i from 1 by 1 to 1001 do > plot_yz := [op(plot_yz), [be_sol2[i,2], be_sol2[i,3]]]: > end do: > plot(plot_yz);

1

0.5

–1

–0.5

0.5 –0.5

–1

148

1


BE, h = 0.001 > be_sol3 := GeneralBackwardEuler(u, f, ic, 0.001, 10): > plot_yz := []: > for i from 1 by 1 to 1001 do > plot_yz := [op(plot_yz), [be_sol3[i,2], be_sol3[i,3]]]: > end do: > plot(plot_yz);

0.6

0.7

0.8

0 –0.2

–0.4

–0.6

–0.8

149

0.9

1


TR, h = 0.1 > tr_sol1 := GeneralTrapezium(u, f, ic, 0.1, 10): > plot_yz := []: > for i from 1 by 1 to 101 do > plot_yz := [op(plot_yz), [tr_sol1[i,2], tr_sol1[i,3]]]: > end do: > plot(plot_yz);

1

0.5

–1

–0.5

0

0.5

–0.5

–1

150

1



1

0.5

–1

–0.5

0.5 –0.5

–1

151

1



0.6

0.7

0.8

0 –0.2

–0.4

–0.6

–0.8

152

0.9

1


B.5.4

Error Assessment

Now compare the analytical solution of this problem to the solutions obtained using the trapezium rule. > ode := diff(diff(y(x),x),x)=-y(x); > ic1 := y(0) = 1; ic2 := D(y)(0) = 0; > trueSoln := dsolve({ode, ic1, ic2}); ode :=

d2 dx2

y(x) = −y(x)

ic1 := y(0) = 1 ic2 := D(y)(0) = 0 trueSoln := y(x) = cos(x) plot_true := []: plot_fe := []: plot_be := []: plot_tr := []: for i from 1 by 1 to 101 do plot_true := [op(plot_true), [fe_sol1[i,1], cos(fe_sol1[i,1])]]: plot_fe := [op(plot_fe), [fe_sol1[i,1], fe_sol1[i,2]]]: plot_be := [op(plot_be), [be_sol1[i,1], be_sol1[i,2]]]: plot_tr := [op(plot_tr), [tr_sol1[i,1], tr_sol1[i,2]]]: end do: plot([plot_true, plot_fe, plot_be, plot_tr], labels=["y", "y’"], legend=["True solution", "Forward Euler’s Method", "Backward Euler’s > Method", "Trapezium Rule Method"]); > > > > > > > > > > > >

1 y’

0.5

2

4

y

6

8

10

0 –0.5 –1 –1.5

True solution Forward Euler’s Method Backward Euler’s Method

153


The plots for the real values and the values obtained using the Trapezium Rule are extremely close, as shown by the graph, which contains two lines, although this is not obvious. > err_plot_fe := []: > err_plot_be := []: > err_plot_tr := []: > for i from 1 by 1 to 101 do > err_plot_fe := [op(err_plot_fe), [fe_sol1[i,1], > abs(cos(fe_sol1[i,1]) - fe_sol1[i,2])]]: > err_plot_be := [op(err_plot_be), [be_sol1[i,1], > abs(cos(be_sol1[i,1]) - be_sol1[i,2])]]: > err_plot_tr := [op(err_plot_tr), [tr_sol1[i,1], > abs(cos(tr_sol1[i,1]) - tr_sol1[i,2])]]: > end do: > plot([err_plot_fe, err_plot_be, err_plot_tr], labels=["x", "Error"], > legend=["Forward Euler’s Method", "Backward Euler’s Method", > "Trapezium Rule Method"]);

0.6 0.5 0.4 Error 0.3 0.2 0.1 0

2

4

x

6

8

10

Forward Euler’s Method Backward Euler’s Method Trapezium Rule Method

For the Forward and Backward Euler Method, the error oscillates with x, producing increasing maxima as x increases. Notably it decreases to 0 for each of the methods where the graph of y crosses the x-axis.

154


We shall now consider the error for each method at x = 2 π (or the closest point returned by the procedure). > evalf(2 * Pi); > fe_err1 := [0.1, abs(fe_sol1[64, 2] - cos(fe_sol1[64, 1]))]; > fe_err2 := [0.01, abs(fe_sol2[629, 2] - cos(fe_sol2[629, 1]))]; > fe_err3 := [0.001, abs(fe_sol3[6284, 2] - cos(fe_sol3[6284, 1]))]; > fe_err := [fe_err1, fe_err2, fe_err3]; > be_err1 := [0.1, abs(be_sol1[64, 2] - cos(be_sol1[64, 1]))]; > be_err2 := [0.01, abs(be_sol2[629, 2] - cos(be_sol2[629, 1]))]; > be_err3 := [0.001, abs(be_sol3[6284, 2] - cos(be_sol3[6284, 1]))]; > be_err := [be_err1, be_err2, be_err3]; > tr_err1 := [0.1, abs(tr_sol1[64, 2] - cos(tr_sol1[64, 1]))]; > tr_err2 := [0.01, abs(tr_sol2[629, 2] - cos(tr_sol2[629, 1]))]; > tr_err3 := [0.001, abs(tr_sol3[6284, 2] - cos(tr_sol3[6284, 1]))]; > tr_err := [tr_err1, tr_err2, tr_err3]; 6.283185308 fe err1 := [0.1, 0.3682471926] fe err2 := [0.01, 0.0318956901] fe err3 := [0.001, 0.0031464292] fe err := [[0.1, 0.3682471926], [0.01, 0.0318956901], [0.001, 0.0031464292]] be err1 := [0.1, 0.2689330814] be err2 := [0.01, 0.0309111282] be err3 := [0.001, 0.0031365719] be err := [[0.1, 0.2689330814], [0.01, 0.0309111282], [0.001, 0.0031365719]] tr err1 := [0.1, 0.0000744028] tr err2 := [0.01, 0.1665 10−6 ] tr err3 := [0.001, 0.5 10−9 ] tr err := [[0.1, 0.0000744028], [0.01, 0.1665 10−6 ], [0.001, 0.5 10−9 ]] We shall now plot the error at x = 2 π against h: > plot([fe_err, be_err, tr_err], legend=["Forward Euler’s Method", > "Backward Euler’s Method", "Trapezium Rule Method"]);

155


0.35 0.3 0.25 0.2 0.15 0.1 0.05 0

0.02

0.04

0.06

0.08

0.1

Forward Euler’s Method Backward Euler’s Method Trapezium Rule Method

The graph demonstrates that the error is proportional to h, as each plot is a straight line.

B.5.5

Conclusion

As h gets smaller the Forward Euler’s Method remains periodic for some time, but eventually spirals outwards, whereas the Backward Euler’s Method spirals inwards. Trapezium Rule Method remains periodic for all h. This is due to the symmetricity property – substituting h → −h into the method leaves it unaltered, so the result for −h is the same as the result for h. This is not the case for the Forward and Backward Euler’s method. Thus the Forward and Backward Euler’s methods are not symmetric. The error is proportional to h. The error oscillates with x, producing increasing maxima as x increases.

156


B.6

4thorder.mws

Fourth Order Rules S. Amen, P.A. Bilokon, A.J. Brinley Codd, M. Fofaria, T. Shah B.6.1

Initialization


B.6.2

Implementation

Implementation : Fourth Order Rules Firstly we implement a basic fourth order method (fourth order symmetric mono implicit Runge-Kutta method) FourthOrderFirst := proc(vars , f, uhalves, fhalves, init, h, interval) local makeeqn, ans, solveEqns, assignList, diffValue, newICSet, i, j, k, ic, icx, ynD, ynDD; ic := init; # get number of variable in u (typically u = [x,y,z]) makeeqn := nops(u); ans := []; # create first x point icx := op(1,ic); ans := [op(ans),ic]; solveEqns := []; 157


assignList := []; # keep looping till get to specified number of # interval points while icx < interval do # icxhalf := icx + h/2; # increment first x by h icx := icx + h; # generate x coordinate solveEqns := [op(1,u) = icx]; #print(solveEqns); assignList := []; # give variables with their initial condition # (or their value at the last iteration) # eg. [x = 1, y = 0, z = 0] for j from 1 by 1 to makeeqn do assignList := [op(assignList), op(j,u) = op(j,ic)]; end do; # evaluate each equation in f # do not evaluate x (because we’ve already incremented that) # using the values in the assign list for i from 2 by 1 to makeeqn - 1 do # get value for yn ynD := eval(op(i,f), assignList); # get value for yn ynDD := eval(op(i+1,f), assignList); # unlike simpler methods, we have two eqns to solve # to find yn1 # list eqns A, B, C (as listed on handout, page 2) solveEqns := [op(solveEqns), op(i,u) = op(i,ic) + h*(ynD) + (h^2/12)*(4*op(i+1,fhalves) + 2*ynDD), op(i, uhalves) = 1/2*(op(i,u) + op(i,ic)) (h^2/48)*(op(i+1,f) + 4*op(i+1,fhalves) + ynDD), 158


op(i+1, u) = ynD + (h/6)*(op(i+1,f) + 4*op(i+1,fhalves) + ynDD)]; end do; # solve for eg. [y,z] and [yhalf,zhalf] etc. newICSet := solve(convert(solveEqns,set), convert([op(u), op(uhalves)],set)); # replace kth "initial condition" with updated value for k from 1 by 1 to makeeqn do ic := subsop(k = eval(op(k,u),newICSet),ic); end do; ans := [op(ans),ic]; end do; # eg. if u = [x,y,z] # matrix returned giving values at each return matrix(ans); end proc: Use Lobatto estimation for half step. FourthOrderFirstLobatto := proc(vars , f, uhalves, fhalves, init, h, interval) local makeeqn, ans, solveEqns, assignList, diffValue, newICSet, i, j, k, ic, icx, ynD, ynDD; ic := init; # get number of variable in u (typically u = [x,y,z]) makeeqn := nops(u); ans := []; # create first x point icx := op(1,ic); ans := [op(ans),ic]; solveEqns := []; assignList := []; # keep looping till get to specified number of 159


# interval points while icx < interval do # icxhalf := icx + h/2; # increment first x by h icx := icx + h; # generate x coordinate solveEqns := [op(1,u) = icx]; assignList := []; # give variables with their initial condition # (or their value at the last iteration) # eg. [x = 1, y = 0, z = 0] for j from 1 by 1 to makeeqn do assignList := [op(assignList), op(j,u) = op(j,ic)]; end do; # evaluate each equation in f # do not evaluate x (because we’ve already incremented that) # using the values in the assign list for i from 2 by 1 to makeeqn - 1 do # get value for yn ynD := eval(op(i,f), assignList); # get value for yn ynDD := eval(op(i+1,f), assignList); # unlike simpler methods, we have two eqns to solve # to find yn1 # gives more accurate estimate for yhalf solveEqns := [ op(solveEqns), op(i,u) = op(i,ic) + h*(ynD) + (h^2/6)*(2*op(i+1,fhalves) + ynDD), op(i, uhalves) = op(i,ic) + (1/2)*h*(ynD) + (h^2/48)* (-op(i+1,f) + 4*op(i+1,fhalves) + 3*ynDD), op(i+1, u) = ynD + (h/6)*(op(i+1,f) + 4*op(i+1,fhalves) + ynDD)]; end do; 160


# find yn1 and ynhalf, by iterative process # solve for eg. [y,z] and [yhalf,zhalf] etc. newICSet := solve(convert(solveEqns,set), convert([op(u),op(uhalves)],set)); # replace kth "initial condition" with updated value for k from 1 by 1 to makeeqn do ic := subsop(k = eval(op(k,u),newICSet),ic); end do; ans := [op(ans),ic]; end do; # eg. if u = [x,y,z] # matrix returned giving values at each return matrix(ans); end proc: Use B approximation for half-step FourthOrderFirstAcc := proc(vars , f, uhalves, fhalves, init, h, interval) local makeeqn, ans, solveEqns, assignList, diffValue, newICSet, i, j, k, ic, icx, ynD, ynDD; ic := init; # get number of variable in u (typically u = [x,y,z]) makeeqn := nops(u); ans := []; # create first x point icx := op(1,ic); ans := [op(ans),ic]; solveEqns := []; assignList := []; # keep looping till get to specified number of # interval points 161


while icx < interval do # increment first x by h icx := icx + h; # generate x coordinate solveEqns := [op(1,u) = icx]; assignList := []; # give variables with their initial condition # (or their value at the last iteration) # eg. [x = 1, y = 0, z = 0] for j from 1 by 1 to makeeqn do assignList := [op(assignList), op(j,u) = op(j,ic)]; end do; # evaluate each equation in f # do not evaluate x (because we’ve already incremented that) # using the values in the assign list for i from 2 by 1 to makeeqn - 1 do # get value for yn ynD := eval(op(i,f), assignList); # get value for yn ynDD := eval(op(i+1,f), assignList); # unlike simpler methods, we have two eqns to solve # to find yn1 # gives more accurate estimate for yhalf eqns (A, B, Z) on # handout solveEqns := [ op(solveEqns), op(i,u) = op(i,ic) + h*(ynD) + (h^2/12)*(4*op(i+1,fhalves) + 2*ynDD), op(i, uhalves) = op(i,ic) +(1/2)*h*(ynD) + (h^2/192)*(-2*op(i+1,f) + 12*op(i+1,fhalves) + 14*ynDD), op(i+1, u) = ynD + (h/6)*(op(i+1,f) + 4*op(i+1,fhalves) + ynDD)]; end do;

162


# solve for eg. [y,z] and [yhalf,zhalf] etc. newICSet := solve(convert(solveEqns,set), convert([op(u),op(uhalves)],set)); # replace kth "initial condition" with updated value for k from 1 by 1 to makeeqn do ic := subsop(k = eval(op(k,u),newICSet),ic); end do; ans := [op(ans),ic]; end do; # eg. if u = [x,y,z] # matrix returned giving values at each return matrix(ans); end proc: Now we can call the above procedures to solve the “circle problem”. Use method A. > u := [x,y,z]; > f := [1,z,-1*y]; > ic := [0,1,0]; > # make a copy of variables that will be used to store ynhalf, > znhalf... > # make copy of f for the same purpose > uh := [xhalf, yhalf, zhalf]; > fh := [1,zhalf,-1*yhalf]; > soln1 := FourthOrderFirst(u,f,uh, fh, ic,0.1,10); u := [x, y, z] f := [1, z, −y] ic := [0, 1, 0] uh := [xhalf , yhalf , zhalf ] fh := [1, zhalf , −yhalf ] 0 0.1000000000 0.2000000000 0.3000000000 0.4000000000 0.5000000000 0.6000000000 0.7000000000 0.8000000000

1 0.9950041666 0.9800665833 0.9553365013 0.9210610155 0.8775825950 0.8253356617 0.7648422496 0.6967067887 163

0 -0.09983340284 -0.1986693036 -0.2955201669 -0.3894182912 -0.4794254777 -0.5646424047 -0.6442176130 -0.7173560136


0.9000000000 1. 1.100000000 1.200000000 1.300000000 1.400000000 1.500000000 1.600000000 1.700000000 1.800000000 1.900000000 2. 2.100000000 2.200000000 2.300000000 2.400000000 2.500000000 2.600000000 2.700000000 2.800000000 2.900000000 3. 3.100000000 3.200000000 3.300000000 3.400000000 3.500000000 3.600000000 3.700000000 3.800000000 3.900000000 4. 4.100000000 4.200000000 4.300000000 4.400000000 4.500000000 4.600000000 4.700000000 4.800000000 4.900000000

0.6216100658 0.5403024223 0.4535962572 0.3623579094 0.2674990022 0.1699673341 0.07073740906 -0.02919930059 -0.1288442605 -0.2272018516 -0.3232893175 -0.4161465844 -0.5048458533 -0.5885008707 -0.6662757835 -0.7373934909 -0.8011434082 -0.8568885677 -0.9040719822 -0.9422222109 -0.9709580692 -0.9899924382 -0.9991351329 -0.9982948023 -0.9874798427 -0.9667983137 -0.9364568582 -0.8967586380 -0.8481003043 -0.7909680351 -0.7259326768 -0.6536440413 -0.5748244124 -0.4902613295 -0.4007997188 -0.3073334509 -0.2107964097 -0.1121531610 -0.01238931529 0.08749832034 0.1865117019

164

-0.7833268320 -0.8414709098 -0.8912072908 -0.9320390256 -0.9635581371 -0.9854496969 -0.9974949718 -0.9995736094 -0.9916648407 -0.9738476875 -0.9463001728 -0.9092975422 -0.8632095136 -0.8084965833 -0.7457054247 -0.6754634260 -0.5984724219 -0.5155016809 -0.4273802189 -0.3349885161 -0.2392497198 -0.1411204200 -0.04158109197 0.05837370046 0.1577452423 0.2555406463 0.3507827733 0.4425199958 0.5298357060 0.6118574745 0.6877657670 0.7568021333 0.8182767849 0.8715754876 0.9161656986 0.9516018872 0.9775299870 0.9936909330 0.9999232505 0.9961646680 0.9824527403


5. 5.100000000 5.200000000 5.300000000 5.400000000 5.500000000 5.600000000 5.700000000 5.800000000 5.900000000 6. 6.100000000 6.200000000 6.300000000 6.400000000 6.500000000 6.600000000 6.700000000 6.800000000 6.900000000 7. 7.100000000 7.200000000 7.300000000 7.400000000 7.500000000 7.600000000 7.700000000 7.800000000 7.900000000 8. 8.100000000 8.200000000 8.300000000 8.400000000 8.500000000 8.600000000 8.700000000 8.800000000 8.900000000 9.

0.2836615207 0.3779770881 0.4685160345 0.5543737247 0.6346922976 0.7086692366 0.7755653888 0.8347123501 0.8855191439 0.9274781256 0.9601700549 0.9832682852 0.9965420266 0.9998586529 0.9931850240 0.9765878214 0.9502328788 0.9143835260 0.8693979578 0.8157256550 0.7539028933 0.6845473853 0.6083521080 0.5260783792 0.4385482506 0.3466362940 0.2512608631 0.1533749174 0.05395650073 -0.04600103131 -0.1454989364 -0.2435430646 -0.3391537917 -0.4313758071 -0.5192876592 -0.6020109622 -0.6787191723 -0.7486458466 -0.8110923012 -0.8654345918 -0.9111297485 165

0.9589244722 0.9258149505 0.8834549944 0.8322678505 0.7727649636 0.7055408668 0.6312672409 0.5506862031 0.4646028923 0.3738774243 0.2794162976 0.1821633365 0.08309026006 -0.01681302666 -0.1165483233 -0.2151191079 -0.3115404941 -0.4048490715 -0.4941125319 -0.5784389846 -0.6569858678 -0.7289683671 -0.7936672575 -0.8504360892 -0.8987076470 -0.9379996175 -0.9679194085 -0.9881680715 -0.9985432884 -0.9989413936 -0.9893584093 -0.9698900856 -0.9407309435 -0.9021723314 -0.8545995140 -0.7984878232 -0.7343979082 -0.6629701341 -0.5849181834 -0.5010219252 -0.4121196229


9.100000000 9.200000000 9.300000000 9.400000000 9.500000000 9.600000000 9.700000000 9.800000000 9.900000000 10.

-0.9477212004 -0.9748433381 -0.9922251660 -0.9996930101 -0.9971722556 -0.9846880882 -0.9623652458 -0.9304267705 -0.8891917811 -0.8390722838

-0.3190995587 -0.2228911581 -0.1244557033 -0.02477672868 0.07514980672 0.1743254703 0.2717593319 0.3664778649 0.4575346731 0.5440199474

Use method B. > u := [x,y,z]; > f := [1,z,-1*y]; > ic := [0,1,0]; > # make a copy of variables that will be used to store ynhalf, > znhalf... > # make copy of f for the same purpose > uh := [xhalf, yhalf, zhalf]; > fh := [1,zhalf,-1*yhalf]; > solnA1 := FourthOrderFirstAcc(u,f,uh, fh, ic,0.1,10); u := [x, y, z] f := [1, z, −y] ic := [0, 1, 0] uh := [xhalf , yhalf , zhalf ] fh := [1, zhalf , −yhalf ] 0 0.1000000000 0.2000000000 0.3000000000 0.4000000000 0.5000000000 0.6000000000 0.7000000000 0.8000000000 0.9000000000 1. 1.100000000 1.200000000 1.300000000 1.400000000

1 0.9950041658 0.9800665799 0.9553364937 0.9210610021 0.8775825743 0.8253356324 0.7648422106 0.6967067390 0.6216100048 0.5403023494 0.4535961722 0.3623578124 0.2674988935 0.1699672144 166

0 -0.09983342015 -0.1986693379 -0.2955202175 -0.3894183571 -0.4794255576 -0.5646424968 -0.6442177155 -0.7173561244 -0.7833269488 -0.8414710301 -0.8912074118 -0.9320391446 -0.9635582513 -0.9854498035


1.500000000 1.600000000 1.700000000 1.800000000 1.900000000 2. 2.100000000 2.200000000 2.300000000 2.400000000 2.500000000 2.600000000 2.700000000 2.800000000 2.900000000 3. 3.100000000 3.200000000 3.300000000 3.400000000 3.500000000 3.600000000 3.700000000 3.800000000 3.900000000 4. 4.100000000 4.200000000 4.300000000 4.400000000 4.500000000 4.600000000 4.700000000 4.800000000 4.900000000 5. 5.100000000 5.200000000 5.300000000 5.400000000 5.500000000

0.07073727921 -0.02919943942 -0.1288444070 -0.2272020039 -0.3232894738 -0.4161467424 -0.5048460107 -0.5885010252 -0.6662759326 -0.7373936318 -0.8011435383 -0.8568886842 -0.9040720825 -0.9422222924 -0.9709581295 -0.9899924750 -0.9991351440 -0.9982947865 -0.9874797979 -0.9667982386 -0.9364567519 -0.8967584999 -0.8481001343 -0.7909678334 -0.7259324442 -0.6536437787 -0.5748241213 -0.4902610119 -0.4007993770 -0.3073330876 -0.2107960280 -0.1121527644 -0.01238890752 0.08749873521 0.1865121196 0.2836619367 0.3779774978 0.4685164331 0.5543741076 0.6346926598 0.7086695734

167

-0.9974950680 -0.9995736926 -0.9916649083 -0.9738477370 -0.9463002020 -0.9092975492 -0.8632094969 -0.8084965415 -0.7457053567 -0.6754633312 -0.5984723001 -0.5155015323 -0.4273800441 -0.3349883162 -0.2392494961 -0.1411201744 -0.04158082669 0.05837398283 0.1577455389 0.2555409539 0.3507830884 0.4425203146 0.5298360246 0.6118577887 0.6877660726 0.7568024261 0.8182770607 0.8715757422 0.9161659279 0.9516020874 0.9775301544 0.9936910642 0.9999233424 0.9961647180 0.9824527462 0.9589244324 0.9258148637 0.8834548598 0.8322676679 0.7727647334 0.7055405899


5.600000000 5.700000000 5.800000000 5.900000000 6. 6.100000000 6.200000000 6.300000000 6.400000000 6.500000000 6.600000000 6.700000000 6.800000000 6.900000000 7. 7.100000000 7.200000000 7.300000000 7.400000000 7.500000000 7.600000000 7.700000000 7.800000000 7.900000000 8. 8.100000000 8.200000000 8.300000000 8.400000000 8.500000000 8.600000000 8.700000000 8.800000000 8.900000000 9. 9.100000000 9.200000000 9.300000000 9.400000000 9.500000000 9.600000000

0.7755656956 0.8347126226 0.8855193778 0.9274783171 0.9601702005 0.9832683816 0.9965420714 0.9998586430 0.9931849586 0.9765876994 0.9502326997 0.9143832900 0.8693976657 0.8157253082 0.7539024939 0.6845469358 0.6083516117 0.5260778401 0.4385476732 0.3466356834 0.2512602248 0.1533742573 0.05395582507 -0.04600171588 -0.1454996229 -0.2435437459 -0.3391544607 -0.4313764566 -0.5192882819 -0.6020115509 -0.6787197201 -0.7486463470 -0.8110927478 -0.8654349788 -0.9111300705 -0.9477214526 -0.9748435161 -0.9922252666 -0.9996930310 -0.9971721948 -0.9846879442

168

0.6312669188 0.5506858380 0.4646024869 0.3738769818 0.2794158218 0.1821628316 0.08308973079 -0.01681357504 -0.1165488852 -0.2151196775 -0.3115410654 -0.4048496383 -0.4941130878 -0.5784395232 -0.6569863827 -0.7289688521 -0.7936677064 -0.8504364961 -0.8987080064 -0.9379999243 -0.9679196580 -0.9881682594 -0.9985434113 -0.9989414486 -0.9893583941 -0.9698899986 -0.9407307839 -0.9021720991 -0.8545992098 -0.7984874486 -0.7343974657 -0.6629696268 -0.5849176152 -0.5010213007 -0.4121189475 -0.3190988384 -0.2228903996 -0.1244549138 -0.02477591574 0.07515063503 0.1743263056


9.700000000 9.800000000 9.900000000 10.

-0.9623650181 -0.9304264600 -0.8891913892 -0.8390718129

0.2717601655 0.3664786880 0.4575354770 0.5440207232

Use Lobatto approximation for the half-step. > u := [x,y,z]; > f := [1,z,-1*y]; > ic := [0,1,0]; > # make a copy of variables that will be used to store ynhalf, > znhalf... > # make copy of f for the same purpose > uh := [xhalf, yhalf, zhalf]; > fh := [1,zhalf,-1*yhalf]; > solnL1 := FourthOrderFirstLobatto(u,f,uh, fh, ic,0.1,10); u := [x, y, z] f := [1, z, −y] ic := [0, 1, 0] uh := [xhalf , yhalf , zhalf ] fh := [1, zhalf , −yhalf ] 0 0.1000000000 0.2000000000 0.3000000000 0.4000000000 0.5000000000 0.6000000000 0.7000000000 0.8000000000 0.9000000000 1. 1.100000000 1.200000000 1.300000000 1.400000000 1.500000000 1.600000000 1.700000000 1.800000000 1.900000000 2.

1 0.9950041666 0.9800665833 0.9553365013 0.9210610155 0.8775825950 0.8253356617 0.7648422496 0.6967067887 0.6216100658 0.5403024223 0.4535962572 0.3623579094 0.2674990022 0.1699673341 0.07073740906 -0.02919930059 -0.1288442605 -0.2272018516 -0.3232893175 -0.4161465844 169

0 -0.09983340284 -0.1986693036 -0.2955201669 -0.3894182912 -0.4794254777 -0.5646424047 -0.6442176130 -0.7173560136 -0.7833268320 -0.8414709098 -0.8912072908 -0.9320390256 -0.9635581371 -0.9854496969 -0.9974949718 -0.9995736094 -0.9916648407 -0.9738476875 -0.9463001728 -0.9092975422


2.100000000 2.200000000 2.300000000 2.400000000 2.500000000 2.600000000 2.700000000 2.800000000 2.900000000 3. 3.100000000 3.200000000 3.300000000 3.400000000 3.500000000 3.600000000 3.700000000 3.800000000 3.900000000 4. 4.100000000 4.200000000 4.300000000 4.400000000 4.500000000 4.600000000 4.700000000 4.800000000 4.900000000 5. 5.100000000 5.200000000 5.300000000 5.400000000 5.500000000 5.600000000 5.700000000 5.800000000 5.900000000 6. 6.100000000

-0.5048458533 -0.5885008707 -0.6662757835 -0.7373934909 -0.8011434082 -0.8568885677 -0.9040719822 -0.9422222109 -0.9709580692 -0.9899924382 -0.9991351329 -0.9982948023 -0.9874798426 -0.9667983136 -0.9364568581 -0.8967586379 -0.8481003042 -0.7909680350 -0.7259326767 -0.6536440412 -0.5748244123 -0.4902613294 -0.4007997187 -0.3073334508 -0.2107964096 -0.1121531609 -0.01238931521 0.08749832040 0.1865117020 0.2836615208 0.3779770882 0.4685160345 0.5543737247 0.6346922976 0.7086692366 0.7755653888 0.8347123501 0.8855191439 0.9274781256 0.9601700549 0.9832682852

170

-0.8632095136 -0.8084965832 -0.7457054246 -0.6754634259 -0.5984724218 -0.5155016808 -0.4273802188 -0.3349885160 -0.2392497197 -0.1411204199 -0.04158109188 0.05837370055 0.1577452424 0.2555406464 0.3507827734 0.4425199959 0.5298357061 0.6118574746 0.6877657671 0.7568021333 0.8182767848 0.8715754875 0.9161656985 0.9516018871 0.9775299869 0.9936909329 0.9999232503 0.9961646680 0.9824527403 0.9589244722 0.9258149505 0.8834549944 0.8322678505 0.7727649636 0.7055408668 0.6312672409 0.5506862031 0.4646028923 0.3738774243 0.2794162976 0.1821633365


6.200000000 6.300000000 6.400000000 6.500000000 6.600000000 6.700000000 6.800000000 6.900000000 7. 7.100000000 7.200000000 7.300000000 7.400000000 7.500000000 7.600000000 7.700000000 7.800000000 7.900000000 8. 8.100000000 8.200000000 8.300000000 8.400000000 8.500000000 8.600000000 8.700000000 8.800000000 8.900000000 9. 9.100000000 9.200000000 9.300000000 9.400000000 9.500000000 9.600000000 9.700000000 9.800000000 9.900000000 10.

0.9965420266 0.9998586529 0.9931850240 0.9765878214 0.9502328788 0.9143835260 0.8693979579 0.8157256551 0.7539028934 0.6845473854 0.6083521081 0.5260783794 0.4385482508 0.3466362942 0.2512608633 0.1533749176 0.05395650092 -0.04600103113 -0.1454989362 -0.2435430645 -0.3391537916 -0.4313758071 -0.5192876593 -0.6020109623 -0.6787191724 -0.7486458467 -0.8110923013 -0.8654345920 -0.9111297487 -0.9477212006 -0.9748433383 -0.9922251662 -0.9996930111 -0.9971722566 -0.9846880892 -0.9623652468 -0.9304267715 -0.8891917821 -0.8390722847

0.08309026007 -0.01681302664 -0.1165483232 -0.2151191078 -0.3115404940 -0.4048490714 -0.4941125318 -0.5784389845 -0.6569858677 -0.7289683670 -0.7936672574 -0.8504360891 -0.8987076469 -0.9379996175 -0.9679194086 -0.9881680716 -0.9985432885 -0.9989413937 -0.9893584095 -0.9698900858 -0.9407309437 -0.9021723316 -0.8545995142 -0.7984878234 -0.7343979084 -0.6629701343 -0.5849181836 -0.5010219254 -0.4121196231 -0.3190995589 -0.2228911583 -0.1244557035 -0.02477672882 0.07514980673 0.1743254704 0.2717593321 0.3664778652 0.4575346736 0.5440199480

We can make a phase plot, for methods A and B for h = 0.1. We see that we do not have any “spiral” shape, like with backward and forward Eulers’ methods. 171

APPENDIX B. MAPLE WORKSHEETS plot_yz := []: plot_yz1 := []: plot_yzA1 := []: for i from 1 by 1 to 101 do # plot_yz := [op(plot_yz),[soln1[i,2], -sin(soln1 [i,1])]]: plot_yz1 := [op(plot_yz1),[soln1[i,2],soln1[i,3]]]: plot_yzA1 := [op(plot_yzA1),[solnA1[i,2],solnA1[i,3]]]: end do: plotsetup("ps", plotoutput="h:/fourth-h01-phase.ps", plotoptions="color"); plot([plot_yz1, plot_yzA1], labels=["y", "y’"], legend=["Fourth Order > A","Fourth Order B"], colour=[red, black], style=[line, line]); > plotsetup(default); > plot([plot_yz1, plot_yzA1], legend=["Fourth Order A","Fourth Order > B"], colour=[red, black], style=[point, line]); > > > > > > > > > > >

1 0.5 –1

–0.5

0.5

1

–0.5 –1 Fourth Order A Fourth Order B

Now compare the analytical solution of this problem to the solutions obtained using the trapezium rule. > ode := diff(diff(y(x),x),x)=-y(x); > ic1 := y(0) = 1; ic2 := D(y)(0) = 0; > trueSoln := dsolve({ode, ic1, ic2});

172


ode :=

d2 dx2

y(x) = −y(x)

ic1 := y(0) = 1 ic2 := D(y)(0) = 0 >

trueSoln := y(x) = cos(x) ode := diff(y(x),‘$‘(x,2)) = -y(x);

>

d ode := dx 2 y(x) = −y(x) trueSoln := y(x) = cos(x);

2

trueSoln := y(x) = cos(x) plot_real := []: plot_estimate := []: plot_accestimate := []: for i from 1 by 1 to 101 do plot_real := [op(plot_real),[soln1[i,1], cos(soln1 [i,1])]]: plot_estimate := [op(plot_estimate),[soln1[i,1], soln1[i,2]]]: plot_accestimate := [op(plot_accestimate),[solnA1[i,1], solnA1[i,2]]]: end do: # plotsetup("ps", plotoutput="h:/fourth-plot.eps", plotoptions="color"); plotsetup(default); plot([plot_estimate, plot_accestimate], legend=["Fourth Order A", "Fourth Order B"] , color=[red,black, green], style=[line,point]); > > > > > > > > > > > > > >

The plots for the real values and the values obtained using the fourth order rule A and fourth order (acc) B are incredibly close, as shown by the graph, which contains two lines, although this is not obvious.

173


1 0.5 0

2

4

6

8

10

–0.5 –1 Fourth Order A Fourth Order B

Error comparison between A and B methods. > err_plot_f := []: > err_plot_fa := []: > err_plot_fl := []: > for i from 1 by 1 to 101 do > err_plot_f := [op(err_plot_f), [soln1[i,1], abs(cos(soln1[i,1]) > soln1[i,2])]]: > err_plot_fa := [op(err_plot_fa), [solnA1[i,1], > abs(cos(solnA1[i,1]) - solnA1[i,2])]]: > err_plot_fl := [op(err_plot_fl), [solnL1[i,1], > abs(cos(solnL1[i,1]) - solnL1[i,2])]]: > end do: > # plotsetup("ps", plotoutput="h:/fourth-errorplot.eps", > plotoptions="color"); > # plot([err_plot_f, err_plot_fa], labels=["x", "Err"], color=[green, > red], style=[line, point], legend=["Fourth Order A", "Fourth Order > B"]); > plotsetup(default); > plot([err_plot_f, err_plot_fa], labels=["x", "Err"], color=[green, > red], style=[line, point], legend=["Fourth Order A", "Fourth Order > B"]); 174


1e–06 8e–07 Err 6e–07 4e–07 2e–07 0

2

4

6

x

8

10


Error comparison between Lobatto and B. > err_plot_f := []: > err_plot_fa := []: > err_plot_fl := []: > for i from 1 by 1 to 101 do > err_plot_f := [op(err_plot_f), [soln1[i,1], abs(cos(soln1[i,1]) > soln1[i,2])]]: > err_plot_fa := [op(err_plot_fa), [solnA1[i,1], > abs(cos(solnA1[i,1]) - solnA1[i,2])]]: > err_plot_fl := [op(err_plot_fl), [solnL1[i,1], > abs(cos(solnL1[i,1]) - solnL1[i,2])]]: > end do: > # plotsetup("ps", plotoutput="h:/fourth-error2plot.eps", > plotoptions="color"); > # plot([err_plot_fa, err_plot_fl], labels=["x", "Err"], color=[green, > red], style=[line, point], legend=["Fourth Order B", "Fourth Order > Lobatto"]); > plotsetup(default); > plot([err_plot_fa, err_plot_fl], labels=["x", "Err"], color=[green, > red], style=[line, point], legend=["Fourth Order B", "Fourth Order > Lobatto"]); 175


1e–06 8e–07 Err 6e–07 4e–07 2e–07 0

2

4

6

x

8

10

Fourth Order B Fourth Order Lobatto

176


B.7

4thorderKepler.mws

Fourth Order Rules For Solving Kepler’s Equation S. Amen, P.A. Bilokon, A.J. Brinley Codd, M. Fofaria, T. Shah B.7.1

Initialization


B.7.2

Implementation

Implementation : Solving Kepler’s Equation using two Fourth Order Rules We have modified the various rules, as we wish to solve Kepler’s equation, which consists of two simultaneous second order differential equations. Essentially it takes in pairs of equations, creating equations for each pair, and then moving on to the next pair. First implement rule A. FourthOrderFirstKepler := proc(vars , f, uhalves, fhalves, svars, init, h, interval) local makeeqn, ans, solveEqns, assignList, diffValue, newICSet, i, j, k, ic, icx, ynD, ynDD, so_far; ic := init; # get number of variable in u (typically u = [x,y,z]) makeeqn := nops(u); ans := []; # create first x point 177


icx := op(1,ic); ans := [op(ans),ic]; solveEqns := []; assignList := []; # keep looping till get to specified number of # interval points while icx < interval do # not really needed # icxhalf := icx + h/2; # increment first x by h icx := icx + h; # generate x coordinate solveEqns := [op(1,u) = icx]; assignList := []; # give variables with their initial condition # (or their value at the last iteration) # eg. [x = 1, y = 0, z = 0] for j from 1 by 1 to makeeqn do assignList := [op(assignList), op(j,u) = op(j,ic)]; end do; # evaluate each equation in f # do not evaluate x (because we’ve already incremented that) # using the values in the assign list # the idea is to examine each pair of equations in turn, adding # them to a list and then solving simultaneously for i from 2 by 2 to makeeqn do # get value for yn differentiated once ynD := eval(op(i,f), assignList); # get value for yn differentiated twice ynDD := eval(op(i+1,f), assignList); # unlike simpler methods, we have two eqns to solve # to find yn1 178


# list eqns A, B, C (as listed on handout, page 2) solveEqns := [ op(solveEqns), op(i,u) = op(i,ic) + h*(ynD) + (h^2/12)*(4*op(i+1,fhalves) + 2*ynDD), op(i, uhalves) = 1/2*(op(i,u) + op(i,ic)) (h^2/48)*(op(i+1,f) + 4*op(i+1,fhalves) + ynDD), op(i+1, u) = ynD + (h/6)*(op(i+1,f) + 4*op(i+1,fhalves) + ynDD)]; end do; # find yn1 and ynhalf, numerically, ie. using fsolve # using solve, to find "exact" solutions takes too long # print(solveEqns); # solve for eg. [y,z] and [yhalf,zhalf] etc. newICSet := fsolve(convert(solveEqns,set), convert([op(u),op(svars)],set)); # replace kth "initial condition" with updated value for k from 1 by 1 to makeeqn do ic := subsop(k = eval(op(k,u),newICSet),ic); end do; ans := [op(ans),ic]; end do; # eg. if u = [x,y,z] # matrix returned giving values at each return matrix(ans); end proc: We can modify the above to get a better estimation for “yhalf”, by using a fourth order method, method B. FourthOrderFirstKeplerAcc := proc(vars , f, uhalves, fhalves, svars, init, h, interval) local makeeqn, ans, solveEqns, assignList, diffValue, newICSet, i, j, k, ic, icx, ynD, ynDD; ic := init;

179


# get number of variable in u (typically u = [x,y,z]) makeeqn := nops(u); ans := []; # create first x point icx := op(1,ic); ans := [op(ans),ic]; solveEqns := []; assignList := []; # keep looping till get to specified number of # interval points while icx < interval do # increment first x by h icx := icx + h; # generate x coordinate solveEqns := [op(1,u) = icx]; # print(solveEqns); assignList := []; # give variables with their initial condition # (or their value at the last iteration) # eg. [x = 1, y = 0, z = 0] for j from 1 by 1 to makeeqn do assignList := [op(assignList), op(j,u) = op(j,ic)]; end do; # evaluate each equation in f # do not evaluate x (because we’ve already incremented that) # using the values in the assign list # the idea is to examine each pair of equations in turn, adding # them to a list and then solving simultaneously for i from 2 by 2 to makeeqn - 1 do # get value for yn differentiated once ynD := eval(op(i,f), assignList); # get value for yn differentiated twice ynDD := eval(op(i+1,f), assignList); 180


# unlike simpler methods, we have two eqns to solve # to find yn1 # gives more accurate estimate for yhalf eqns A, B, Z solveEqns := [ op(solveEqns), op(i,u) = op(i,ic) + h*(ynD) + (h^2/12)*(4*op(i+1,fhalves) + 2*ynDD), op(i, uhalves) = op(i,ic) +(1/2)*h*(ynD) + (h^2/192)*(-2*op(i+1,f) + 12*op(i+1,fhalves) + 14*ynDD), op(i+1, u) = ynD + (h/6)*(op(i+1,f) + 4*op(i+1,fhalves) + ynDD)]; end do; # solve for eg. [y,z] and [yhalf,zhalf] etc. newICSet := fsolve(convert(solveEqns,set), convert([op(u),op(svars)],set)); # replace kth "initial condition" with updated value for k from 1 by 1 to makeeqn do ic := subsop(k = eval(op(k,u),newICSet),ic); end do; ans := [op(ans),ic]; end do; # eg. if u = [x,y,z] # matrix returned giving values at each return matrix(ans); end proc: Let us begin by describing Kepler’s Equation. We consider y(x) and z(x) to be components of direction as functions of x which is essentially time. Thus solving the ODE and getting y(x) and z(x) will enable us to plot the orbit of an object as time varies. y(x) d2 y(x) = − 3 2 2 dx (y(x) + z(x)2 )( 2 ) z(x) d2 z(x) = − 3 2 dx (y(x)2 + z(x)2 )( 2 ) 181

APPENDIX B. MAPLE WORKSHEETS with the initial conditions (where we choose e = 0.6) y(0) = 1 − e, d z(0) = 1+e z(0) = 0, dx 1−e .

d dx y(0)

= 0,

We can use the same trick that we used earlier to make into a set of four first order differential equations. d y(x) = a(x) dx y(x) d a(x) = − 3 dx (y(x)2 + z(x)2 )( 2 ) d z(x) = b(x) dx z(x) d b(x) = − 3 2 dx (y(x) + z(x)2 )( 2 ) with the initial conditions being modified to take into account new vars (where we 1+e choose e = 0.6) y(0) = 1 − e, a(0) = 0, z(0) = 0, b(0) = 1−e . We can then adjust h as too increase accuracy. Try to solve analytically. Maple does not produce a solution! We could solve using polar coordinates though. > ode1 := diff(diff(y(x),x),x)=-y(x)/(y(x)^2+z(x)^2)^(3/2): > ode2 := diff(diff(z(x),x),x)=-z(x)/(y(x)^2+z(x)^2)^(3/2): > ecc := 0.6: > ic1 := y(0) = 1-ecc; ic2 := D(y)(0) = 0; ic3 := z(0) = 0; ic4 := > D(z)(0) = sqrt((1+ecc)/(1-ecc)): > trueSoln := dsolve({ode1, ode2, ic1, ic2, ic3, ic4}, {y(x), > z(x)}): Solve Kepler’s Equation using consts (problem described on http://bartok.ucsc.edu/peter/115/math intro/node22.html) - this problem find the rotation of the earth around the sun

182

APPENDIX B. MAPLE WORKSHEETS u := [x,y,a,z,b]: ecc := 0.5: GM := 4*Pi^2: f := [1,a,-GM*y/((y^2+z^2)^(3/2)),b,-GM*z/((y^2+z^2)^(3/2))]: # initial conditions ic := [0,0,2*Pi,1,-Pi]: # make a copy of variables that will be used to store ynhalf, znhalf... # make copy of f for the same purpose uh := [xhalf, yhalf, ahalf, zhalf, bhalf]: sv := [yhalf, zhalf]: fh := [1,ahalf,-GM*yhalf/((yhalf^2+zhalf^2)^(3/2)),bhalf,-GM*zhalf/((yhalf^2 +zhalf^2)^(3/2))]: solnKsA1 := FourthOrderFirstKeplerAcc(u,f,uh, fh, sv, ic, 0.1, 10): > solnKs1 := FourthOrderFirstKepler(u,f,uh, fh, sv, ic, 0.1, 10): > > > > > > > > > > > > > > >

We now attempt to solve Kepler’s Equations as described in our project. > u := [x,y,a,z,b]: > ecc := 0.6: > f := [1,a,-y/((y^2+z^2)^(3/2)),b,-z/((y^2+z^2)^(3/2))]: > # initial conditions > ic := [0,1-ecc,0,0,((1+ecc)/(1-ecc))^(1/2)]: > # make a copy of variables that will be used to store ynhalf, > znhalf... > # make copy of f for the same purpose > uh := [xhalf, yhalf, ahalf, zhalf, bhalf]: > sv := [yhalf, zhalf]: > fh := > [1,ahalf,-yhalf/((yhalf^2+zhalf^2)^(3/2)),bhalf,-zhalf/((yhalf^2+zhalf > ^2)^(3/2))]: > solnKA1 := FourthOrderFirstKeplerAcc(u, f, uh, fh, sv, ic, 0.1, 10); > solnK1 := FourthOrderFirstKepler(u, f, uh, fh, sv, ic, 0.1, 10); 0 0.1000000000 0.2000000000 0.3000000000 0.4000000000 0.5000000000 0.6000000000 0.7000000000 0.8000000000 0.9000000000 1.000000000

0.4 0.3698493407 0.2899170806 0.1808150702 0.05905311296 -0.06575964731 -0.1886747505 -0.3072750181 -0.4204410193 -0.5277136429 -0.6289705762

0 -0.5831814464 -0.9788548839 -1.175527933 -1.244096738 -1.244209938 -1.210226409 -1.159996061 -1.102600756 -1.042669491 -0.9825610602 183

0 0.1950024383 0.3649299360 0.4998220431 0.6016711559 0.6762448393 0.7291563387 0.7649106767 0.7869397947 0.7978383971 0.7995851098

2.000000000 1.855583318 1.527357437 1.175084108 0.8719912489 0.6289302191 0.4367902863 0.2839804549 0.1608952773 0.06035133687 -0.02288245909


1.100000000 1.200000000 1.300000000 1.400000000 1.500000000 1.600000000 1.700000000 1.800000000 1.900000000 2.000000000 2.100000000 2.200000000 2.300000000 2.400000000 2.500000000 2.600000000 2.700000000 2.800000000 2.900000000 3.000000000 3.100000000 3.200000000 3.300000000 3.400000000 3.500000000 3.600000000 3.700000000 3.800000000 3.900000000 4.000000000 4.100000000 4.200000000 4.300000000 4.400000000 4.500000000 4.600000000 4.700000000 4.800000000 4.900000000 5.000000000 5.100000000 5.200000000

-0.7242601287 -0.8137142387 -0.8975021684 -0.9758056283 -1.048805449 -1.116674591 -1.179574679 -1.237654504 -1.291049613 -1.339882498 -1.384263079 -1.424289328 -1.460047931 -1.491614934 -1.519056351 -1.542428698 -1.561779471 -1.577147551 -1.588563530 -1.596049975 -1.599621620 -1.599285495 -1.595040983 -1.586879817 -1.574786008 -1.558735709 -1.538697008 -1.514629661 -1.486484748 -1.454204254 -1.417720591 -1.376956044 -1.331822162 -1.282219104 -1.228034983 -1.169145261 -1.105412298 -1.036685252 -0.9628006278 -0.8835840371 -0.7988541308 -0.7084304190

-0.9234559601 -0.8659116129 -0.8101515844 -0.7562199169 -0.7040652181 -0.6535871918 -0.6046626059 -0.5571597723 -0.5109465038 -0.4658943228 -0.4218805048 -0.3787888686 -0.3365098495 -0.2949401690 -0.2539822878 -0.2135437498 -0.1735364804 -0.1338760735 -0.09448108563 -0.05527234372 -0.01617226856 0.02289578935 0.06200820516 0.1012417329 0.1406741616 0.1803849812 0.2204560609 0.2609723434 0.3020225521 0.3436999022 0.3861027967 0.4293354673 0.4735084920 0.5187390705 0.5651508531 0.6128729799 0.6620377456 0.7127758869 0.7652077574 0.8194273327 0.8754735590 0.9332789815 184

0.7937105870 0.7814180161 0.7636683423 0.7412407701 0.7147761846 0.6848087461 0.6517892203 0.6161024601 0.5780806945 0.5380137709 0.4961571555 0.4527382639 0.4079615364 0.3620125584 0.3150614510 0.2672657005 0.2187725546 0.1697210850 0.1202439970 0.07046924666 0.02052152113 -0.02947637440 -0.07940218229 -0.1291328732 -0.1785433058 -0.2275048265 -0.2758837701 -0.3235398136 -0.3703241295 -0.4160772732 -0.4606267207 -0.5037839563 -0.5453409719 -0.5850660038 -0.6226982711 -0.6579413977 -0.6904550813 -0.7198444054 -0.7456459398 -0.7673094128 -0.7841731938 -0.7954310077

-0.09261430626 -0.1516437881 -0.2020569405 -0.2454306779 -0.2829739568 -0.3156259319 -0.3441249856 -0.3690579103 -0.3908954224 -0.4100181431 -0.4267358383 -0.4413018272 -0.4539238811 -0.4647725349 -0.4739874669 -0.4816824154 -0.4879489683 -0.4928594715 -0.4964692334 -0.4988181534 -0.4999318665 -0.4998224639 -0.4984888305 -0.4959166134 -0.4920778230 -0.4869300435 -0.4804152129 -0.4724579045 -0.4629630146 -0.4518127203 -0.4388625204 -0.4239360999 -0.4068186640 -0.3872482444 -0.3649042873 -0.3393925413 -0.3102248427 -0.2767917655 -0.2383251564 -0.1938461350 -0.1420919518 -0.08141177004


5.300000000 5.400000000 5.500000000 5.600000000 5.700000000 5.800000000 5.900000000 6.000000000 6.100000000 6.200000000 6.300000000 6.400000000 6.500000000 6.600000000 6.700000000 6.800000000 6.900000000 7.000000000 7.100000000 7.200000000 7.300000000 7.400000000 7.500000000 7.600000000 7.700000000 7.800000000 7.900000000 8.000000000 8.100000000 8.200000000 8.300000000 8.400000000 8.500000000 8.600000000 8.700000000 8.800000000 8.900000000 9.000000000 9.100000000 9.200000000 9.300000000 9.400000000

-0.6121480961 -0.5098856607 -0.4016163614 -0.2875050639 -0.1680938540 -0.04466437491 0.08004830822 0.2004872892 0.3060051637 0.3789387162 0.3991141531 0.3593199572 0.2730070520 0.1607904559 0.03797196589 -0.08680415104 -0.2091374696 -0.3268949059 -0.4391019118 -0.5453735381 -0.6456254726 -0.7399251811 -0.8284137481 -0.9112641916 -0.9886590614 -1.060778464 -1.127793829 -1.189864861 -1.247138279 -1.299747542 -1.347813103 -1.391442928 -1.430733138 -1.465768667 -1.496623909 -1.523363313 -1.546041914 -1.564705794 -1.579392475 -1.590131236 -1.596943372 -1.599842365

0.9925767208 1.052729223 1.112406130 1.168963350 1.217217663 1.246995242 1.238321797 1.153097086 0.9274363338 0.4949889489 -0.1051503637 -0.6655608476 -1.024495865 -1.194381719 -1.247982017 -1.240479276 -1.202777727 -1.150820622 -1.092715421 -1.032595973 -0.9725653795 -0.9136693995 -0.8563932507 -0.8009214441 -0.7472769480 -0.6953972292 -0.6451764069 -0.5964887867 -0.5492019625 -0.5031839811 -0.4583070929 -0.4144495277 -0.3714961304 -0.3293383449 -0.2878738348 -0.2470059094 -0.2066428548 -0.1666972276 -0.1270851411 -0.08772556017 -0.04853961157 -0.009449909076 185

-0.8000860824 -0.7968871282 -0.7842380419 -0.7600703216 -0.7216660610 -0.6654295014 -0.5866614907 -0.4796102282 -0.3387490197 -0.1634501216 0.03358639239 0.2258026807 0.3900898137 0.5190869976 0.6159026230 0.6864407546 0.7361535127 0.7693540873 0.7893174785 0.7985213806 0.7988592022 0.7917994790 0.7784996515 0.7598865131 0.7367134013 0.7096013328 0.6790690058 0.6455550090 0.6094345025 0.5710319290 0.5306308313 0.4884815384 0.4448072609 0.3998089874 0.3536694684 0.3065565028 0.2586256860 0.2100227445 0.1608855505 0.1113458952 0.06153107956 0.01156537593

-0.009616557537 0.07623836607 0.1801664385 0.3076854219 0.4663572954 0.6661968587 0.9188806032 1.231946793 1.587720066 1.897668124 1.995591099 1.808207055 1.466543275 1.119720770 0.8268163341 0.5930945467 0.4083248327 0.2611091219 0.1422583725 0.04495344486 -0.03576499348 -0.1035113737 -0.1609483192 -0.2100645021 -0.2523668068 -0.2890131062 -0.3209047385 -0.3487516766 -0.3731190642 -0.3944608925 -0.4131446874 -0.4294698281 -0.4436812892 -0.4559800503 -0.4665310444 -0.4754692646 -0.4829044709 -0.4889248169 -0.4935996298 -0.4969815099 -0.4991078726 -0.5000020185


9.500000000 9.600000000 9.700000000 9.800000000 9.900000000 10.00000000

-1.598834009 -1.593916454 -1.585080191 -1.572307975 -1.555574671 -1.534847047

0 0.1000000000 0.2000000000 0.3000000000 0.4000000000 0.5000000000 0.6000000000 0.7000000000 0.8000000000 0.9000000000 1.000000000 1.100000000 1.200000000 1.300000000 1.400000000 1.500000000 1.600000000 1.700000000 1.800000000 1.900000000 2.000000000 2.100000000 2.200000000 2.300000000 2.400000000 2.500000000 2.600000000 2.700000000 2.800000000 2.900000000 3.000000000 3.100000000 3.200000000

0.4 0.3698428971 0.2898942069 0.1807716631 0.05898930502 -0.06584176971 -0.1887725482 -0.3073859472 -0.4205628666 -0.5278445773 -0.6291091298 -0.7244051552 -0.8138648694 -0.8976577725 -0.9759657795 -1.048969898 -1.116843241 -1.179747570 -1.237831799 -1.291231586 -1.340069522 -1.384455622 -1.424487949 -1.460253272 -1.491827723 -1.519277396 -1.542658888 -1.562019779 -1.577399032 -1.588827324 -1.596327310 -1.599913818 -1.599593974

0.02962010894 -0.03842859650 0.06874683410 -0.08832854243 0.1080071415 -0.1380111702 0.1474790471 -0.1873508442 0.1872423745 -0.2362181709 0.2273794361 -0.2844784782 Table B.5: solnKA1

0 -0.5833124717 -0.9790517251 -1.175738236 -1.244292312 -1.244380162 -1.210370035 -1.160115696 -1.102700158 -1.042752454 -0.9826310247 -0.9235159195 -0.8659641363 -0.8101988724 -0.7562638640 -0.7041074699 -0.6536291927 -0.6047056398 -0.5572049956 -0.5109949723 -0.4659470141 -0.4219383366 -0.3788527143 -0.3365805516 -0.2950185508 -0.2540691634 -0.2136399337 -0.1736427962 -0.1339933622 -0.09461021404 -0.05541421242 -0.01632782030 0.02272556072 186

0 0.1950018222 0.3649265445 0.4998122390 0.6016514225 0.6762127663 0.7291107362 0.7648512954 0.7868670136 0.7977529863 0.7994880708 0.7936030495 0.7813011723 0.7635434060 0.7411089500 0.7146386677 0.6846666876 0.6516437368 0.6159546256 0.5779315379 0.5378642743 0.4960082529 0.4525908405 0.4078164272 0.3618705474 0.3149232696 0.2671320258 0.2186440069 0.1695982253 0.1201273230 0.07035918865 0.02041843665 -0.02957220702

-0.4996737872 -0.4981198319 -0.4953235274 -0.4912545072 -0.4858678055 -0.4791025583

2.000000000 1.855569226 1.527311820 1.175001023 0.8718778324 0.6287990827 0.4366524852 0.2838437542 0.1607646399 0.06022976603 -0.02299324058 -0.09271339254 -0.1517307906 -0.2021317977 -0.2454935350 -0.2830250892 -0.3156656971 -0.3441537921 -0.3690761988 -0.3909036538 -0.4100167920 -0.4267253897 -0.4412827760 -0.4538967325 -0.4647378066 -0.4739456929 -0.4816341500 -0.4878947914 -0.4927999950 -0.4964051085 -0.4987500796 -0.4998606015 -0.4997488364


3.300000000 3.400000000 3.500000000 3.600000000 3.700000000 3.800000000 3.900000000 4.000000000 4.100000000 4.200000000 4.300000000 4.400000000 4.500000000 4.600000000 4.700000000 4.800000000 4.900000000 5.000000000 5.100000000 5.200000000 5.300000000 5.400000000 5.500000000 5.600000000 5.700000000 5.800000000 5.900000000 6.000000000 6.100000000 6.200000000 6.300000000 6.400000000 6.500000000 6.600000000 6.700000000 6.800000000 6.900000000 7.000000000 7.100000000 7.200000000 7.300000000 7.400000000

-1.595367264 -1.587225527 -1.575152892 -1.559125636 -1.539111984 -1.515071839 -1.486956437 -1.454707939 -1.418258948 -1.377531956 -1.332438739 -1.282879708 -1.228743249 -1.169905124 -1.106228015 -1.037561424 -0.9637422105 -0.8845963348 -0.7999427543 -0.7096011747 -0.6134067463 -0.5112374382 -0.4030650108 -0.2890509526 -0.1697302183 -0.04636978295 0.07832501786 0.1988553329 0.3046729674 0.3782192283 0.3992589567 0.3602573370 0.2744157688 0.1623774477 0.03956597123 -0.08528516342 -0.2077274666 -0.3256043948 -0.4379304378 -0.5443157816 -0.6446742076 -0.7390727085

0.06182224549 0.1010389176 0.1404532852 0.1801447451 0.2201950603 0.2606890531 0.3017153104 0.3433668936 0.3857420333 0.4289447694 0.4730854683 0.5182810999 0.5646550716 0.6123362775 0.6614567834 0.7121471527 0.7645276910 0.8186925790 0.8746814563 0.9324284856 0.9916701539 1.051775613 1.111427706 1.168008416 1.216386610 1.246493255 1.238560769 1.154853459 0.9318791780 0.5027401064 -0.09621835176 -0.6590997527 -1.021450907 -1.193645548 -1.248430946 -1.241453535 -1.203945818 -1.152024984 -1.093884123 -1.033698652 -0.9735918460 -0.9146191842 187

-0.07949057168 -0.1292137244 -0.1786166313 -0.2275707599 -0.2759425825 -0.3235919334 -0.3703701663 -0.4161180470 -0.4606632984 -0.5038176957 -0.5453735770 -0.5850995941 -0.6227354692 -0.6579854407 -0.6905099654 -0.7199150735 -0.7457385295 -0.7674315889 -0.7843345978 -0.7956438780 -0.8003661190 -0.7972547154 -0.7847199919 -0.7607023547 -0.7224962312 -0.6665224538 -0.5881018424 -0.4814974664 -0.3411598210 -0.1663274240 0.03055772795 0.2230524798 0.3878426702 0.5173326803 0.6145423640 0.6853734912 0.7352990182 0.7686531958 0.7887274014 0.7980111930 0.7984064297 0.7913876840

-0.4984137538 -0.4958411036 -0.4920030192 -0.4868572335 -0.4803458638 -0.4723937009 -0.4629059060 -0.4517649802 -0.4388268208 -0.4239156068 -0.4068171596 -0.3872702865 -0.3649554187 -0.3394795678 -0.3103562075 -0.2769780599 -0.2385798247 -0.1941864639 -0.1425404748 -0.08199829099 -0.01038109390 0.07524128451 0.1788615619 0.3059682427 0.4640841006 0.6631803485 0.9149206298 1.226996009 1.582411168 1.894120406 1.996349026 1.812610620 1.471786622 1.124200149 0.8302230921 0.5955851755 0.4101252139 0.2624080453 0.1431953608 0.04562819878 -0.03528147879 -0.1031681394


7.500000000 7.600000000 7.700000000 7.800000000 7.900000000 8.000000000 8.100000000 8.200000000 8.300000000 8.400000000 8.500000000 8.600000000 8.700000000 8.800000000 8.900000000 9.000000000 9.100000000 9.200000000 9.300000000 9.400000000 9.500000000 9.600000000 9.700000000 9.800000000 9.900000000 10.00000000

-0.8276525791 -0.9105873577 -0.9880602374 -1.060251998 -1.127334727 -1.189468748 -1.246801361 -1.299466558 -1.347585283 -1.391265958 -1.430605124 -1.465688107 -1.496589672 -1.523374614 -1.546098299 -1.564807130 -1.579538938 -1.590323306 -1.597181828 -1.600128287 -1.599168777 -1.594301755 -1.585518028 -1.572800677 -1.556124908 -1.535457844

-0.8572704308 -0.8017320512 -0.7480276691 -0.6960947657 -0.6458271731 -0.5970987852 -0.5497767504 -0.5037286828 -0.4588264342 -0.4149478792 -0.3719775534 -0.3298066368 -0.2883325712 -0.2474584842 -0.2070925170 -0.1671471155 -0.1275383144 -0.08818503032 -0.04900837097 -0.009930959485 0.02912372806 0.06823201726 0.1074706877 0.1469176291 0.1866525065 0.2267574390 Table B.6: solnK1

188

0.7781167541 0.7595236188 0.7363639760 0.7092606126 0.6787335702 0.6452224661 0.6091032552 0.5707009973 0.5302997169 0.4881501191 0.4444757074 0.3994776969 0.3533390109 0.3062275767 0.2582990807 0.2096993077 0.1605661602 0.1110314335 0.06122240843 0.01126331337 -0.03872330030 -0.08861522990 -0.1382893030 -0.1876200324 -0.2364782058 -0.2847293703

-0.1607083399 -0.2099005172 -0.2522584846 -0.2889451017 -0.3208653480 -0.3487318944 -0.3731119112 -0.3944609326 -0.4131476759 -0.4294724519 -0.4436809744 -0.4559748179 -0.4665194026 -0.4754501275 -0.4828770982 -0.4888887701 -0.4935547406 -0.4969278603 -0.4990457840 -0.4999320487 -0.4995967360 -0.4980367547 -0.4952357570 -0.4911636852 -0.4857759249 -0.4790120191


We can plot the yz plane. > # nb: includes commented out subsections to > facilitate creation of EPS files for graphs > plot_real := []: > plot_Kestimate := []: > plot_Kaccestimate := []: > plot_Ksestimate := []: > plot_Ksaccestimate := []: > for i from 1 by 1 to 101 do > # to plot Earth Kepler problem > plot_Ksestimate := [op(plot_Ksestimate),[solnKs1[i,4], > solnKs1[i,2]]]: > plot_Ksaccestimate := [op(plot_Ksaccestimate),[solnKsA1[i,4], > solnKsA1[i,2]]]: > # to plot Kepler Equation (described in project, without equations) > plot_Kestimate := [op(plot_Kestimate),[solnK1[i,4], solnK1[i,2]]]: > plot_Kaccestimate := [op(plot_Kaccestimate),[solnKA1[i,4], > solnKA1[i,2]]]: > end do: > # plot Earth Kepler > # plotsetup("ps", plotoutput="h:/fourth-kepler-earth001.eps", > plotoptions="color"); > plotsetup(default); > plot([plot_Ksestimate, plot_Ksaccestimate],labels=["z","y"], > legend=["Fourth Order A Kepler Earth Problem", "Fourth Order B Kepler > Earth Problem"], style=[line, line]); > # plot Kepler problem described in notes > # plotsetup("ps", plotoutput="h:/fourth-kepler.eps", > plotoptions="color"); > plotsetup(default); > plot([plot_Kestimate, plot_Kaccestimate], labels=["z","y"], > legend=["Fourth Order A Kepler Problem", "Fourth Order B Kepler > Problem"], style=[line, line]);

189


–1

–0.5

0.5 0.5

z

1

1.5

0 –0.5 y

–1 –1.5 –2 Fourth Order A Kepler Earth Problem Fourth Order B Kepler Earth Problem

–0.8 –0.6 –0.4 –0.2

0.2

z 0.4

0.6

0.8

0 –0.5 y –1 –1.5

Fourth Order A Kepler Problem Fourth Order B Kepler Problem

190


B.8

modtrap.mws

Modified Trapezium Rule (with comparison with Fourth Order A) S. Amen, P.A. Bilokon, A.J. Brinley Codd, M. Fofaria, T. Shah Version 1.0

Version History • 10 October 2004. Version 1.0. Complete Functionality

B.8.1

Initialization


B.8.2

Implementation

Implementation : Modified Trapezium Rule We first implement a modified trapezium rule, then follow with a fourth order rule (A) that we can use for comparison. ModTrap := proc(vars , f, init, h, interval) local makeeqn, ans, solveEqns, assignList, diffValue, newICSet, i, j, k, ic, icx, ynD, ynDD; ic := init; # get number of variable in u (typically u = [x,y,z]) 191


makeeqn := nops(u); ans := []; # create first x point icx := op(1,ic); ans := [op(ans),ic]; solveEqns := []; assignList := []; # keep looping till get to specified number of # interval points while icx < interval do # increment first x by h icx := icx + h; # generate x coordinate solveEqns := [op(1,u) = icx]; assignList := []; # give variables with their initial condition # (or their value at the last iteration) # eg. [x = 1, y = 0, z = 0] for j from 1 by 1 to makeeqn do assignList := [op(assignList), op(j,u) = op(j,ic)]; end do; # evaluate each equation in f # do not evaluate x (because we’ve already incremented that) # using the values in the assign list for i from 2 by 1 to makeeqn - 1 do # get value for yn ynD := eval(op(i,f), assignList); # get value for yn ynDD := eval(op(i+1,f), assignList); # unlike simpler methods, we have two eqns to solve # to find yn1 # now solve 192


solveEqns := [ op(solveEqns), op(i,u) = op(i,ic) + h*(ynD) + (h^2/4)*(op(i+1,f) + ynDD), op(i+1, u) = ynD + (h/2)*(op(i+1,f) + ynDD)]; end do; # find yn1 and ynhalf, by iterative process # solve for eg. [y,z] and [yhalf,zhalf] etc. newICSet := solve(convert(solveEqns,set), convert([op(u)],set)); # replace kth "initial condition" with updated value for k from 1 by 1 to makeeqn do ic := subsop(k = eval(op(k,u),newICSet),ic); end do; ans := [op(ans),ic]; end do; # eg. if u = [x,y,z] # matrix returned giving values at each return matrix(ans); end proc: ModTrapKepler := proc(vars , f, init, h, interval) local makeeqn, ans, solveEqns, assignList, diffValue, newICSet, i, j, k, ic, icx, ynD, ynDD, so_far; ic := init; # get number of variable in u (typically u = [x,y,z]) makeeqn := nops(u); ans := []; # create first x point icx := op(1,ic); ans := [op(ans),ic];

193


solveEqns := []; assignList := []; # keep looping till get to specified number of # interval points while icx < interval do # not really needed # icxhalf := icx + h/2; # increment first x by h icx := icx + h; # generate x coordinate solveEqns := [op(1,u) = icx]; assignList := []; # give variables with their initial condition # (or their value at the last iteration) # eg. [x = 1, y = 0, z = 0] for j from 1 by 1 to makeeqn do assignList := [op(assignList), op(j,u) = op(j,ic)]; end do; # evaluate each equation in f # do not evaluate x (because we’ve already incremented that) # using the values in the assign list # the idea is to examine each pair of equations in turn, adding # them to a list and then solving simultaneously for i from 2 by 2 to makeeqn do # get value for yn differentiated once ynD := eval(op(i,f), assignList); # get value for yn differentiated twice ynDD := eval(op(i+1,f), assignList); # unlike simpler methods, we have two eqns to solve # to find yn1 solveEqns := [ op(solveEqns), 194


op(i,u) = op(i,ic) + h*(ynD) + (h^2/4)*(op(i+1,f) + ynDD), op(i+1, u) = ynD + (h/2)*(op(i+1,f) + ynDD)]; end do; # find yn1 and ynhalf, numerically, ie. using fsolve # using solve, to find "exact" solutions takes too long # solve for eg. [y,z] and [yhalf,zhalf] etc. newICSet := fsolve(convert(solveEqns,set),convert([op(u)],set)); # replace kth "initial condition" with updated value for k from 1 by 1 to makeeqn do ic := subsop(k = eval(op(k,u),newICSet),ic); end do; ans := [op(ans),ic]; end do; # eg. if u = [x,y,z] # matrix returned giving values at each return matrix(ans); end proc: Fourth order A, for Kepler FourthOrderFirstKepler := proc(vars , f, uhalves, fhalves, svars, init, h, interval) local makeeqn, ans, solveEqns, assignList, diffValue, newICSet, i, j, k, ic, icx, ynD, ynDD, so_far; ic := init; # get number of variable in u (typically u = [x,y,z]) makeeqn := nops(u); ans := []; # create first x point icx := op(1,ic); ans := [op(ans),ic]; solveEqns := []; 195


assignList := []; # keep looping till get to specified number of # interval points while icx < interval do # not really needed # icxhalf := icx + h/2; # increment first x by h icx := icx + h; # generate x coordinate solveEqns := [op(1,u) = icx]; assignList := []; # give variables with their initial condition # (or their value at the last iteration) # eg. [x = 1, y = 0, z = 0] for j from 1 by 1 to makeeqn do assignList := [op(assignList), op(j,u) = op(j,ic)]; end do; # evaluate each equation in f # do not evaluate x (because we’ve already incremented that) # using the values in the assign list # the idea is to examine each pair of equations in turn, adding # them to a list and then solving simultaneously for i from 2 by 2 to makeeqn do # get value for yn differentiated once ynD := eval(op(i,f), assignList); # get value for yn differentiated twice ynDD := eval(op(i+1,f), assignList); # unlike simpler methods, we have two eqns to solve # to find yn1 # list eqns A, B, C (as listed on handout, page 2) solveEqns := [ op(solveEqns), 196


op(i,u) = op(i,ic) + h*(ynD) + (h^2/12)*(4*op(i+1,fhalves) + 2*ynDD), op(i, uhalves) = 1/2*(op(i,u) + op(i,ic)) (h^2/48)*(op(i+1,f) + 4*op(i+1,fhalves) + ynDD), op(i+1, u) = ynD + (h/6)*(op(i+1,f) + 4*op(i+1,fhalves) + ynDD)]; end do; # solve for eg. [y,z] and [yhalf,zhalf] etc. newICSet := fsolve(convert(solveEqns,set), convert([op(u), op(svars)],set)); # replace kth "initial condition" with updated value for k from 1 by 1 to makeeqn do ic := subsop(k = eval(op(k,u),newICSet),ic); end do; ans := [op(ans),ic]; end do; # eg. if u = [x,y,z] # matrix returned giving values at each return matrix(ans); end proc: FourthOrderFirst := proc(vars , f, uhalves, fhalves, init, h, interval) local makeeqn, ans, solveEqns, assignList, diffValue, newICSet, i, j, k, ic, icx, ynD, ynDD; ic := init; # get number of variable in u (typically u = [x,y,z]) makeeqn := nops(u); ans := []; # create first x point icx := op(1,ic); ans := [op(ans),ic]; solveEqns := []; 197


assignList := []; # keep looping till get to specified number of # interval points while icx < interval do # icxhalf := icx + h/2; # increment first x by h icx := icx + h; # generate x coordinate solveEqns := [op(1,u) = icx]; #print(solveEqns); assignList := []; # give variables with their initial condition # (or their value at the last iteration) # eg. [x = 1, y = 0, z = 0] for j from 1 by 1 to makeeqn do assignList := [op(assignList), op(j,u) = op(j,ic)]; end do; # evaluate each equation in f # do not evaluate x (because we’ve already incremented that) # using the values in the assign list for i from 2 by 1 to makeeqn - 1 do # get value for yn ynD := eval(op(i,f), assignList); # get value for yn ynDD := eval(op(i+1,f), assignList); # unlike simpler methods, we have two eqns to solve # to find yn1 solveEqns := [ op(solveEqns), op(i,u) = op(i,ic) + h*(ynD) + (h^2/12)*(4*op(i+1,fhalves) + 2*ynDD), op(i, uhalves) = 1/2*(op(i,u) + op(i,ic)) (h^2/48)*(op(i+1,f) + 4*op(i+1,fhalves) + ynDD), op(i+1, u) = ynD + (h/6)*(op(i+1,f) + 198


4*op(i+1,fhalves) + ynDD) ]; end do; # print(solveEqns); # find yn1 and ynhalf, by iterative process # solve for eg. [y,z] and [yhalf,zhalf] etc. newICSet := solve(convert(solveEqns,set), convert([op(u),op(uhalves)],set)); #print(newICSet); # replace kth "initial condition" with updated value for k from 1 by 1 to makeeqn do ic := subsop(k = eval(op(k,u),newICSet),ic); end do; ans := [op(ans),ic]; end do; # eg. if u = [x,y,z] # matrix returned giving values at each return matrix(ans); end proc: Try solving circle problem using fourth order A. > u := [x,y,z]: f := [1,z,-1*y]: ic := [0,1,0]: # make a copy of variables that will be used to store ynhalf, znhalf... # make copy of f for the same purpose uh := [xhalf, yhalf, zhalf]: fh := [1,zhalf,-1*yhalf]: solnF1 := FourthOrderFirst(u,f,uh, fh, ic,0.1,10); 0 0.1000000000 0.2000000000 0.3000000000 0.4000000000 0.5000000000 0.6000000000

1 0.9950041666 0.9800665833 0.9553365013 0.9210610155 0.8775825950 0.8253356617 199

0 -0.09983340284 -0.1986693036 -0.2955201669 -0.3894182912 -0.4794254777 -0.5646424047


0.7000000000 0.8000000000 0.9000000000 1. 1.100000000 1.200000000 1.300000000 1.400000000 1.500000000 1.600000000 1.700000000 1.800000000 1.900000000 2. 2.100000000 2.200000000 2.300000000 2.400000000 2.500000000 2.600000000 2.700000000 2.800000000 2.900000000 3. 3.100000000 3.200000000 3.300000000 3.400000000 3.500000000 3.600000000 3.700000000 3.800000000 3.900000000 4. 4.100000000 4.200000000 4.300000000 4.400000000 4.500000000 4.600000000 4.700000000

0.7648422496 0.6967067887 0.6216100658 0.5403024223 0.4535962572 0.3623579094 0.2674990022 0.1699673341 0.07073740906 -0.02919930059 -0.1288442605 -0.2272018516 -0.3232893175 -0.4161465844 -0.5048458533 -0.5885008707 -0.6662757835 -0.7373934909 -0.8011434082 -0.8568885677 -0.9040719822 -0.9422222109 -0.9709580692 -0.9899924382 -0.9991351329 -0.9982948023 -0.9874798427 -0.9667983137 -0.9364568582 -0.8967586380 -0.8481003043 -0.7909680351 -0.7259326768 -0.6536440413 -0.5748244124 -0.4902613295 -0.4007997188 -0.3073334509 -0.2107964097 -0.1121531610 -0.01238931529

200

-0.6442176130 -0.7173560136 -0.7833268320 -0.8414709098 -0.8912072908 -0.9320390256 -0.9635581371 -0.9854496969 -0.9974949718 -0.9995736094 -0.9916648407 -0.9738476875 -0.9463001728 -0.9092975422 -0.8632095136 -0.8084965833 -0.7457054247 -0.6754634260 -0.5984724219 -0.5155016809 -0.4273802189 -0.3349885161 -0.2392497198 -0.1411204200 -0.04158109197 0.05837370046 0.1577452423 0.2555406463 0.3507827733 0.4425199958 0.5298357060 0.6118574745 0.6877657670 0.7568021333 0.8182767849 0.8715754876 0.9161656986 0.9516018872 0.9775299870 0.9936909330 0.9999232505


4.800000000 4.900000000 5. 5.100000000 5.200000000 5.300000000 5.400000000 5.500000000 5.600000000 5.700000000 5.800000000 5.900000000 6. 6.100000000 6.200000000 6.300000000 6.400000000 6.500000000 6.600000000 6.700000000 6.800000000 6.900000000 7. 7.100000000 7.200000000 7.300000000 7.400000000 7.500000000 7.600000000 7.700000000 7.800000000 7.900000000 8. 8.100000000 8.200000000 8.300000000 8.400000000 8.500000000 8.600000000 8.700000000 8.800000000

0.08749832034 0.1865117019 0.2836615207 0.3779770881 0.4685160345 0.5543737247 0.6346922976 0.7086692366 0.7755653888 0.8347123501 0.8855191439 0.9274781256 0.9601700549 0.9832682852 0.9965420266 0.9998586529 0.9931850240 0.9765878214 0.9502328788 0.9143835260 0.8693979578 0.8157256550 0.7539028933 0.6845473853 0.6083521080 0.5260783792 0.4385482506 0.3466362940 0.2512608631 0.1533749174 0.05395650073 -0.04600103131 -0.1454989364 -0.2435430646 -0.3391537917 -0.4313758071 -0.5192876592 -0.6020109622 -0.6787191723 -0.7486458466 -0.8110923012

201

0.9961646680 0.9824527403 0.9589244722 0.9258149505 0.8834549944 0.8322678505 0.7727649636 0.7055408668 0.6312672409 0.5506862031 0.4646028923 0.3738774243 0.2794162976 0.1821633365 0.08309026006 -0.01681302666 -0.1165483233 -0.2151191079 -0.3115404941 -0.4048490715 -0.4941125319 -0.5784389846 -0.6569858678 -0.7289683671 -0.7936672575 -0.8504360892 -0.8987076470 -0.9379996175 -0.9679194085 -0.9881680715 -0.9985432884 -0.9989413936 -0.9893584093 -0.9698900856 -0.9407309435 -0.9021723314 -0.8545995140 -0.7984878232 -0.7343979082 -0.6629701341 -0.5849181834


8.900000000 9. 9.100000000 9.200000000 9.300000000 9.400000000 9.500000000 9.600000000 9.700000000 9.800000000 9.900000000 10.

-0.8654345918 -0.9111297485 -0.9477212004 -0.9748433381 -0.9922251660 -0.9996930101 -0.9971722556 -0.9846880882 -0.9623652458 -0.9304267705 -0.8891917811 -0.8390722838

-0.5010219252 -0.4121196229 -0.3190995587 -0.2228911581 -0.1244557033 -0.02477672868 0.07514980672 0.1743254703 0.2717593319 0.3664778649 0.4575346731 0.5440199474

Solve using modified trapezium rule. > u := [x,y,z]: > f := [1,z,-1*y]: > ic := [0,1,0]: > solnMT1 := ModTrap(u,f, ic,0.1,10); 0 0.1000000000 0.2000000000 0.3000000000 0.4000000000 0.5000000000 0.6000000000 0.7000000000 0.8000000000 0.9000000000 1. 1.100000000 1.200000000 1.300000000 1.400000000 1.500000000 1.600000000 1.700000000 1.800000000 1.900000000 2. 2.100000000 2.200000000

1 0.9950124688 0.9800996262 0.9554102287 0.9211905546 0.8777819473 0.8256174104 0.7652172884 0.6971840762 0.6221964094 0.5410022947 0.4544116483 0.3632882174 0.2685409639 0.1711149975 0.07198214832 -0.02786872728 -0.1274416106 -0.2257432560 -0.3217930983 -0.4146330344 -0.5033369800 -0.5870201079 202

0 -0.09975062344 -0.1985062282 -0.2952817209 -0.3891117600 -0.4790603851 -0.5642303530 -0.6437720879 -0.7168921561 -0.7828611804 -0.8410211156 -0.8907918127 -0.9316768060 -0.9632682651 -0.9852510632 -0.9974059205 -0.9996115916 -0.9918460747 -0.9741868314 -0.9468100137 -0.9099887071 -0.8640902064 -0.8095723520


2.300000000 2.400000000 2.500000000 2.600000000 2.700000000 2.800000000 2.900000000 3. 3.100000000 3.200000000 3.300000000 3.400000000 3.500000000 3.600000000 3.700000000 3.800000000 3.900000000 4. 4.100000000 4.200000000 4.300000000 4.400000000 4.500000000 4.600000000 4.700000000 4.800000000 4.900000000 5. 5.100000000 5.200000000 5.300000000 5.400000000 5.500000000 5.600000000 5.700000000 5.800000000 5.900000000 6. 6.100000000 6.200000000 6.300000000

-0.6648476736 -0.7360433423 -0.7998969328 -0.8557715015 -0.9031096962 -0.9414393154 -0.9703780187 -0.9896371407 -0.9990245706 -0.9984466683 -0.9879091980 -0.9675172719 -0.9374743006 -0.8980799648 -0.8497272254 -0.7928984039 -0.7281603715 -0.6561588940 -0.5776121905 -0.4933037694 -0.4040746125 -0.3108147861 -0.2144545628 -0.1159551419 -0.01629906126 0.08351960354 0.1825051551 0.2796702064 0.3740455300 0.4646897261 0.5506986132 0.6312142473 0.7054334799 0.7726159696 0.8320915668 0.8832669987 0.9256317872 0.9587633409 0.9823311705 0.9961001856 0.9999330394

203

-0.7469789629 -0.6769344121 -0.6001373984 -0.5173539767 -0.4294099168 -0.3371824662 -0.2415915995 -0.1435908416 -0.04415775603 0.05571580592 0.1550335992 0.2528049227 0.3480545013 0.4398322145 0.5272225740 0.6093538555 0.6854067943 0.7546227576 0.8163113118 0.8698571098 0.9147260289 0.9504704988 0.9767339662 0.9932544514 0.9998671616 0.9965061348 0.9832048968 0.9600961287 0.9274103419 0.8854735791 0.8347041621 0.7756085190 0.7087761326 0.6348736601 0.5546382833 0.4688703551 0.3784254158 0.2842056594 0.1871509339 0.08822936612 -0.01157229513


6.400000000 6.500000000 6.600000000 6.700000000 6.800000000 6.900000000 7. 7.100000000 7.200000000 7.300000000 7.400000000 7.500000000 7.600000000 7.700000000 7.800000000 7.900000000 8. 8.100000000 8.200000000 8.300000000 8.400000000 8.500000000 8.600000000 8.700000000 8.800000000 8.900000000 9. 9.100000000 9.200000000 9.300000000 9.400000000 9.500000000 9.600000000 9.700000000 9.800000000 9.900000000 10.

0.9937914986 0.9777368256 0.9519291668 0.9166259552 0.8721793425 0.8190326865 0.7577161285 0.6888413048 0.6130952462 0.5312335243 0.4440727148 0.3524822524 0.2573757575 0.1597019234 0.06043505263 -0.03943466156 -0.1389110126 -0.2370017176 -0.3327283157 -0.4251359281 -0.5133027831 -0.5963494108 -0.6734474160 -0.7438277412 -0.8067883383 -0.8617011715 -0.9080184818 -0.9452782511 -0.9731088109 -0.9912325496 -0.9994686823 -0.9977350514 -0.9860489519 -0.9645269526 -0.9333837369 -0.8929299602 -0.8435691514

-0.1112585220 -0.2098349382 -0.3063182378 -0.3997459939 -0.4891862588 -0.5737468602 -0.6525843009 -0.7249121725 -0.7900090000 -0.8472254385 -0.8959907504 -0.9358184987 -0.9663113992 -0.9871652833 -0.9981721321 -0.9992221519 -0.9903048682 -0.9715092317 -0.9430227300 -0.9051295178 -0.8582075822 -0.8027249725 -0.7392351312 -0.6683713733 -0.5908405693 -0.5074160938 -0.4189301111 -0.3262652744 -0.2303459213 -0.1321288533 -0.03259379171 0.06726639498 0.1664555952 0.2639843904 0.3588799248 0.4501956096 0.5370205652

We can make a phase plot, for both methods for h =0.1, alongside phase plot from the exact solution. We see that we do not have any “spiral” shape, like with backward and forward Eulers’ methods.

204

APPENDIX B. MAPLE WORKSHEETS plot_yz := []: plot_yz1 := []: plot_yzMT1 := []: for i from 1 by 1 to 101 do #plot_yz := [op(plot_yz),[soln1[i,2], -sin(soln1 [i,1])]]: plot_yz1 := [op(plot_yz1),[solnF1[i,2],solnF1[i,3]]]: plot_yzMT1 := [op(plot_yzMT1),[solnMT1[i,2],solnMT1[i,3]]]: end do: #plotsetup("ps", plotoutput="h:/modtrap-phase.ps", plotoptions="color"); #plot([plot_yz1, plot_yzMT1], labels=["y", "y’"], legend=["Fourth Order A","Mod Trap"], colour=[red, black], style=[line, line]); plotsetup(default); plot([plot_yz1, plot_yzMT1], legend=["Fourth Order A","Mod Trap"], > colour=[red, black], style=[line, line]); > > > > > > > > > > > > > >

1 0.5 –1

–0.5

0.5

1

–0.5 –1 Fourth Order A Mod Trap

Now compare the analytical solution of this problem to the solutions obtained using the trapezium rule. > ode := diff(diff(y(x),x),x)=-y(x): > ic1 := y(0) = 1; ic2 := D(y)(0) = 0: > trueSoln := dsolve({ode, ic1, ic2}); ic1 := y(0) = 1 205


>

trueSoln := y(x) = cos(x) ode := diff(y(x),‘$‘(x,2)) = -y(x);

>

d ode := dx 2 y(x) = −y(x) trueSoln := y(x) = cos(x);

2

trueSoln := y(x) = cos(x) plot_real := []: plot_estimate := []: plot_MTestimate := []: for i from 1 by 1 to 101 do plot_real := [op(plot_real),[solnF1[i,1], cos(solnF1 [i,1])]]: plot_estimate := [op(plot_estimate),[solnF1[i,1], solnF1[i,2]]]: plot_MTestimate := [op(plot_MTestimate),[solnMT1[i,1], solnMT1[i,2]]]: end do: #plotsetup("ps", plotoutput="h:/modtrap-fourth-plot.eps", plotoptions="color"); plotsetup(default); plot([plot_real, plot_estimate, plot_MTestimate], legend=["Real Sol", > "Fourth Order A", "Mod Trap"] , color=[red,black, green], > style=[line,line,point]); > > > > > > > > > > > > >

The plots for the real values and the values obtained using the fourth order rule A and modified trapezium rule are incredibly close, as shown by the graph, which contains two lines, although this is not obvious.

1 0.5 0

2

4

6

8

–0.5 –1 Real Sol Fourth Order A Mod Trap

206

10


Comparing errors between modified trapezium and fourth order A. We see fourth order A produces a far smaller error. > err_plot_f := []: > err_plot_f4 := []: > err_plot_fMT := []: > for i from 1 by 1 to 101 do > #err_plot_f := [op(err_plot_f), [soln1[i,1], abs(cos(soln1[i,1]) > soln1[i,2])]]: > err_plot_f4 := [op(err_plot_f4), [solnF1[i,1], > abs(cos(solnF1[i,1]) - solnF1[i,2])]]: > err_plot_fMT := [op(err_plot_fMT), [solnMT1[i,1], > abs(cos(solnMT1[i,1]) - solnMT1[i,2])]]: > end do: > # plotsetup("ps", plotoutput="mt-errorplot.eps", plotoptions="color"); > # plot([err_plot_f4, err_plot_fMT], labels=["x", "Err"], color=[green, > red], style=[line, point], legend=["Fourth Order A", "Modified > Trapezium"]); > plotsetup(default); > plot([err_plot_f4, err_plot_fMT], labels=["x", "Err"], color=[green, > red], style=[point, point], legend=["Fourth Order A", "Modified > Trapezium"]);

0.006 0.005 0.004 Err

0.003 0.002 0.001 0

2

4

x

6

8

Fourth Order A Modified Trapezium

207

10


Solve Kepler’s Equation problem as described in notes using Modified Trapezium Rule. > u := [x,y,a,z,b]; > ecc := 0.6; > f := [1,a,-y/((y^2+z^2)^(3/2)),b,-z/((y^2+z^2)^(3/2))]; > ic := [0,1-ecc,0,0,((1+ecc)/(1-ecc))^(1/2)]: > ic := [0,1-ecc,0,0,((1+ecc)/(1-ecc))^(1/2)]: > solnMTK1 := ModTrapKepler(u, f, ic, 0.1, 10); u := [x, y, a, z, b] ecc := 0.6 y z f := [1, a, − 2 , b, − 2 ] (3/2) 2 (y + z ) (y + z 2 )(3/2) 0 , 0.4 , 0 , 0 , 2.000000000 0.1000000000 , 0.3717461794 , -0.5650764113 , 0.1934289064 , 1.868578127 0.2000000000 , 0.2954734956 , -0.9603772650 , 0.3649025651 , 1.560895046 0.3000000000 , 0.1892688510 , -1.163715628 , 0.5041069212 , 1.223192076 0.4000000000 , 0.06912615463 , -1.239138300 , 0.6118008007 , 0.9306855150 0.5000000000 , -0.05511826243 , -1.245750041 , 0.6931618727 , 0.6965359268 0.6000000000 , -0.1783235145 , -1.218355001 , 0.7536046312 , 0.5123192433 0.7000000000 , -0.2979920172 , -1.175015052 , 0.7975492753 , 0.3665736390 0.8000000000 , -0.4129844678 , -1.124833959 , 0.8283644842 , 0.2497305393 0.9000000000 , -0.5228458659 , -1.072394002 , 0.8485836330 , 0.1546524370 1.000000000 , -0.6274654111 , -1.019996902 , 0.8601246637 , 0.07616817656 1.100000000 , -0.7269040261 , -0.9687753985 , 0.8644597702 , 0.01053395236 1.200000000 , -0.8213054183 , -0.9192524448 , 0.8627370569 , -0.04498821867 1.300000000 , -0.9108494994 , -0.8716291768 , 0.8558661873 , -0.09242917303 1.400000000 , -0.9957277982 , -0.8259367988 , 0.8445788149 , -0.1333182744 1.500000000 , -1.076130555 , -0.7821183497 , 0.8294716716 , -0.1688245928 1.600000000 , -1.152240141 , -0.7400733803 , 0.8110377154 , -0.1998545312 1.700000000 , -1.224227937 , -0.6996825433 , 0.7896889817 , -0.2271201431 1.800000000 , -1.292253118 , -0.6608210863 , 0.7657735946 , -0.2511876002 1.900000000 , -1.356462479 , -0.6233661195 , 0.7395886113 , -0.2725120662 2.000000000 , -1.416990803 , -0.5872003563 , 0.7113898509 , -0.2914631405 2.100000000 , -1.473961513 , -0.5522138484 , 0.6813995102 , -0.3083436711 2.200000000 , -1.527487435 , -0.5183045857 , 0.6498121346 , -0.3234038398 2.300000000 , -1.577671588 , -0.4853784655 , 0.6167993509 , -0.3368518341 2.400000000 , -1.624607957 , -0.4533489245 , 0.5825136582 , -0.3488620218 2.500000000 , -1.668382222 , -0.4221364040 , 0.5470914933 , -0.3595812771 2.600000000 , -1.709072429 , -0.3916677479 , 0.5106557334 , -0.3691339215 2.700000000 , -1.746749596 , -0.3618755896 , 0.4733177564 , -0.3776256179 208


2.800000000 , -1.781478263 , -0.3326977574 , 2.900000000 , -1.813316987 , -0.3040767166 , 3.000000000 , -1.842318776 , -0.2759590529 , 3.100000000 , -1.868531478 , -0.2482949994 , 3.200000000 , -1.891998128 , -0.2210380065 , 3.300000000 , -1.912757246 , -0.1941443511 , 3.400000000 , -1.930843102 , -0.1675727815 , 3.500000000 , -1.946285950 , -0.1412841947 , 3.600000000 , -1.959112226 , -0.1152413413 , 3.700000000 , -1.969344721 , -0.08940855527 , 3.800000000 , -1.977002725 , -0.06375150446 , 3.900000000 , -1.982102148 , -0.03823695827 , 4.000000000 , -1.984655625 , -0.01283257028 , 4.100000000 , -1.984672588 , 0.01249332751 , 4.200000000 , -1.982159326 , 0.03777192148 , 4.300000000 , -1.977119025 , 0.06303410643 , 4.400000000 , -1.969551786 , 0.08831067138 , 4.500000000 , -1.959454629 , 0.1136324831 , 4.600000000 , -1.946821472 , 0.1390306697 , 4.700000000 , -1.931643098 , 0.1645368064 , 4.800000000 , -1.913907102 , 0.1901831051 , 4.900000000 , -1.893597816 , 0.2160026098 , 5.000000000 , -1.870696216 , 0.2420294002 , 5.100000000 , -1.845179806 , 0.2682988043 , 5.200000000 , -1.817022485 , 0.2948476222 , 5.300000000 , -1.786194386 , 0.3217143620 , 5.400000000 , -1.752661694 , 0.3489394882 , 5.500000000 , -1.716386436 , 0.3765656823 , 5.600000000 , -1.677326247 , 0.4046381138 , 5.700000000 , -1.635434106 , 0.4332047170 , 5.800000000 , -1.590658046 , 0.4623164659 , 5.900000000 , -1.542940840 , 0.4920276328 , 6.000000000 , -1.492219658 , 0.5223960076 , 6.100000000 , -1.438425705 , 0.5534830406 , 6.200000000 , -1.381483860 , 0.5853538482 , 6.300000000 , -1.321312318 , 0.6180769849 , 6.400000000 , -1.257822278 , 0.6517238278 , 6.500000000 , -1.190917720 , 0.6863673254 , 6.600000000 , -1.120495368 , 0.7220797079 , 6.700000000 , -1.046444958 , 0.7589284999 , 6.800000000 , -0.9686500460 , 0.7969697388 , 6.900000000 , -0.8869897317 , 0.8362365464 , 209

0.4351791524 , -0.3851464624 0.3963331565 , -0.3917734571 0.3568658587 , -0.3975724997 0.3168572340 , -0.4025999934 0.2763820266 , -0.4069041549 0.2355105148 , -0.4105260799 0.1943091801 , -0.4135006129 0.1528412966 , -0.4158570568 0.1111674562 , -0.4176197502 0.06934604201 , -0.4188085337 0.02743365921 , -0.4194391224 -0.01451446677 , -0.4195233973 -0.05644411792 , -0.4190696256 -0.09830173008 , -0.4180826177 -0.1400340521 , -0.4165638243 -0.1815878122 , -0.4145113797 -0.2229093857 , -0.4119200902 -0.2639444586 , -0.4087813687 -0.3046376827 , -0.4050831135 -0.3449323147 , -0.4008095279 -0.3847698348 , -0.3959408750 -0.4240895365 , -0.3904531599 -0.4628280810 , -0.3843177295 -0.5009190064 , -0.3775007763 -0.5382921817 , -0.3699627293 -0.5748731937 , -0.3616575098 -0.6105826505 , -0.3525316242 -0.6453353846 , -0.3425230576 -0.6790395336 , -0.3315599219 -0.7115954696 , -0.3195587987 -0.7428945445 , -0.3064226989 -0.7728176063 , -0.2920385369 -0.8012332324 , -0.2762739850 -0.8279956081 , -0.2589735280 -0.8529419584 , -0.2399534766 -0.8758894129 , -0.2189956118 -0.8966311440 , -0.1958390101 -0.9149315657 , -0.1701694246 -0.9305203042 , -0.1416053445 -0.9430845457 , -0.1096794870 -0.9522592165 , -0.07381393017 -0.9576142271 , -0.03328628290


7.000000000 , -0.8013419112 , 0.8767198650 , -0.9586376948 , 0.01281692970 7.100000000 , -0.7115891305 , 0.9183357498 , -0.9547135828 , 0.06566531028 7.200000000 , -0.6176288871 , 0.9608691181 , -0.9450914973 , 0.1267764005 7.300000000 , -0.5193916668 , 1.003875288 , -0.9288453572 , 0.1981464003 7.400000000 , -0.4168727101 , 1.046503846 , -0.9048162037 , 0.2824366702 7.500000000 , -0.3101887377 , 1.087175603 , -0.8715325483 , 0.3832364390 7.600000000 , -0.1996812468 , 1.122974216 , -0.8270999739 , 0.5054150491 7.700000000 , -0.08610890714 , 1.148472577 , -0.7690530193 , 0.6555240432 7.800000000 , 0.02898679944 , 1.153441554 , -0.6941774736 , 0.8419868716 7.900000000 , 0.1425834185 , 1.118490829 , -0.5983772598 , 1.074017405 8.000000000 , 0.2489036508 , 1.007913817 , -0.4768896105 , 1.355735579 8.100000000 , 0.3375229254 , 0.7644716751 , -0.3257672192 , 1.666712246 8.200000000 , 0.3924386206 , 0.3338422294 , -0.1462290487 , 1.924051164 8.300000000 , 0.3969599624 , -0.2434153920 , 0.04922997932 , 1.985129396 8.400000000 , 0.3454005357 , -0.7877731407 , 0.2377580465 , 1.785431948 8.500000000 , 0.2489165615 , -1.141906343 , 0.3986385918 , 1.432178957 8.600000000 , 0.1267080668 , -1.302263553 , 0.5238914962 , 1.072879129 8.700000000 , -0.005482678215 , -1.341551348 , 0.6162326017 , 0.7739429804 8.800000000 , -0.1385515180 , -1.319825448 , 0.6819855344 , 0.5411156742 8.900000000 , -0.2680664600 , -1.270473393 , 0.7271326902 , 0.3618274417 9.000000000 , -0.3920886763 , -1.209970934 , 0.7563517665 , 0.2225540842 9.100000000 , -0.5098957360 , -1.146170259 , 0.7731159655 , 0.1127298981 9.200000000 , -0.6213415986 , -1.082746992 , 0.7799894189 , 0.02473917173 9.300000000 , -0.7265468540 , -1.021358117 , 0.7788856697 , -0.04681415628 9.400000000 , -0.8257485462 , -0.9626757282 , 0.7712562944 , -0.1057733507 9.500000000 , -0.9192267953 , -0.9068892543 , 0.7582219770 , -0.1549129966 9.600000000 , -1.007268979 , -0.8539544133 , 0.7406628752 , -0.1962690397 9.700000000 , -1.090152639 , -0.8037188059 , 0.7192813973 , -0.2313605180 9.800000000 , -1.168137875 , -0.7559859053 , 0.6946464168 , -0.2613390913 9.900000000 , -1.241464550 , -0.7105475806 , 0.6672249104 , -0.2870910364 10.00000000 , -1.310351942 , -0.6672002575 , 0.6374049669 , -0.3093078341 Solve Kepler’s Equation using Fourth Order A.

210

APPENDIX B. MAPLE WORKSHEETS > > > > > > > > > > > > > >

u := [x,y,a,z,b]; ecc := 0.6; f := [1,a,-y/((y^2+z^2)^(3/2)),b,-z/((y^2+z^2)^(3/2))]; ic := [0,1-ecc,0,0,((1+ecc)/(1-ecc))^(1/2)]; # make a copy of variables that will be used to store ynhalf, znhalf... # make copy of f for the same purpose # sv is the vars we are interested in plotting uh := [xhalf, yhalf, ahalf, zhalf, bhalf]; sv := [yhalf, zhalf]; fh := [1,ahalf,-yhalf/((yhalf^2+zhalf^2)^(3/2)),bhalf,-zhalf/((yhalf^2+zhalf ^2)^(3/2))]; solnK1 := FourthOrderFirstKepler(u, f, uh, fh, sv, ic, 0.1, 10): u := [x, y, a, z, b] ecc := 0.6 z y , b, − 2 ] f := [1, a, − 2 (3/2) 2 (y + z ) (y + z 2 )(3/2) ic := [0, 0.4, 0, 0, 2.000000000] uh := [xhalf , yhalf , ahalf , zhalf , bhalf ] sv := [yhalf , zhalf ] fh := [1, ahalf , −

zhalf yhalf ] 2 (3/2) , bhalf , − 2 (yhalf + zhalf ) (yhalf + zhalf 2 )(3/2) 2

Plot Kepler’s Equation result from modified Trapezium and from fourth order A. > plot_real := []: > plot_Kestimate := []: > plot_Ktrapestimate := []: > for i from 1 by 1 to 101 do > # plot_real := [op(plot_real),[soln[i,1], cos(soln [i,1])]]: > plot_Kestimate := [op(plot_Kestimate),[solnK1[i,4], solnK1[i,2]]]: > plot_Ktrapestimate := [op(plot_Ktrapestimate),[solnMTK1[i,4], > solnMTK1[i,2]]]: > end do: > # plot second problem as described in notes > # plotsetup("ps", plotoutput="h:/fourth-kepler.eps", > plotoptions="color"); > plotsetup(default); > plot([plot_Kestimate, plot_Ktrapestimate], labels=["z","y"], > legend=["Fourth Order A Kepler Problem", "Modified Trapezium Kepler > Problem"], style=[line, line]);

211


z 0.2 0.4 0.6 0.8

–0.8 –0.6 –0.4 –0.2 0

–0.5 –1

y

–1.5 –2 Fourth Order A Kepler Problem Modified Trapezium Kepler Problem

212


B.9

stiff.mws

Solving Stiff Equations S. Amen, P.A. Bilokon, A.J. Brinley Codd, M. Fofaria, T. Shah B.9.1

Initialization

The following commands initialize Maple. The “warnings” that appear have no adverse effects on the execution, so can be safely ignored. > restart; > with(DEtools): > with(plots): > with(linalg): Warning, the name changecoords has been redefined Warning, the name adjoint has been redefined Warning, the protected names norm and trace have been redefined and unprotected

B.9.2

Background

Introduction The equation we wish to solve a second order ODE

(

d d2 y(x)) + (λ + 1) ( y(x)) + λ y = 0 2 dx dx

with the initial conditions y(0) = 1,

d dx

y(0) = −1.

By introducing auxillary variables we can convert it into a set of first order ODEs, where λ is a constant. We shall investigate how varying λ affects the numerical solution. . d y(x) = z(x) dx d z(x) = −λ y(x) − (λ + 1) z(x) dx

213


The above is an example of a stiff equation. Essentially this is where solutions can change rapidly with the respect to the independent variable (here x), yet the other solution changes a lot more slowly. Thus a very small. If we solve the above analytically, we find that: y(x) = e(−x) We also have an unwanted solution where y(x) = e(−λ x) Issues We aim to demonstrate that under certain circumstances forward Euler fails to converge to the real solution, and instead blows up to infinity. In fact forward Euler is only stable for the step sizes that satisfy h
u := [x,y,z]; > f := [1,z,-la103*y - (la103 + 1)*z]; > ic := [0,1,-1]; u := [x, y, z] f := [1, z, −1000 y − 1001 z] ic := [0, 1, −1] >

la103 := 10^3; la103 := 1000

If we limit the range for x = [0, 1.5], there doesn’t seem to be a problem, until x = 1.5 where we have a sudden jump in value.

214

APPENDIX B. MAPLE WORKSHEETS soln1h01x15 := GeneralBackwardEuler(u,f,ic,0.1,1.5): soln2h01x15 := GeneralForwardEuler(u,f,ic,0.1,1.5): plot_yz := []: plot_yzLa := []: plot_yz1 := []: plot_yz2 := []: for i from 1 by 1 to 16 do plot_yz := [op(plot_yz),[soln1h01x15[i,1], exp(-soln1h01x15 [i,1])]]: plot_yzLa := [op(plot_yzLa),[soln1h01x15[i,1], exp(-la103*(soln1h01x15 [i,1]))]]: plot_yz1 := [op(plot_yz1),[soln1h01x15[i,1], soln1h01x15[i,2]]]: plot_yz2 := [op(plot_yz2),[soln2h01x15[i,1],soln2h01x15[i,2]]]: end do: plotsetup("ps", plotoutput="h:\stiff-ploth01.eps", plotoptions="color"); plot([plot_yz, plot_yzLa, plot_yz1, plot_yz2], legend=["Real sol", > "Unwanted sol", "Backward Euler","Forward Euler"], colour=[green, > brown, red, black], style=[line,line,line, point]); > plotsetup(default); > plot([plot_yz, plot_yzLa, plot_yz1, plot_yz2], legend=["Real sol", > "Unwanted sol", "Backward Euler","Forward Euler"], colour=[green, > brown, red, black], style=[line,line,line, point]); > > > > > > > > > > > > > > > > >

1 0.8 0.6 0.4 0.2 0

0.2

0.4

0.6

0.8

1

1.2

Real sol Unwanted sol Backward Euler

215

1.4


However, when increasing the range to x = 1.6, we find that forward Euler breaks down far more significantly. Increasing x further, and forward Euler becomes very unstable, blowing up to infinity. > soln1h01x16 := GeneralBackwardEuler(u,f,ic,0.1,1.6): > soln2h01x16 := GeneralForwardEuler(u,f,ic,0.1,1.6): > plot_yz := []: > plot_yzLa := []: > plot_yz1 := []: > plot_yz2 := []: > for i from 1 by 1 to 17 do > plot_yz := [op(plot_yz),[soln1h01x16[i,1], exp(-soln1h01x16 > [i,1])]]: > plot_yzLa := [op(plot_yzLa),[soln1h01x16[i,1], exp(-la*(soln1h01x16 > [i,1]))]]: > plot_yz1 := [op(plot_yz1),[soln1h01x16[i,1], soln1h01x16[i,2]]]: > plot_yz2 := [op(plot_yz2),[soln2h01x16[i,1],soln2h01x16[i,2]]]: > end do: > plotsetup("ps", plotoutput="h:\stiff-ploth0116.eps", > plotoptions="color"); > plot([plot_yz, plot_yzLa, plot_yz1, plot_yz2], legend=["Real sol", > "Unwanted sol", "Backward Euler","Forward Euler"], colour=[green, > brown, red, black], style=[line,line,line, point]); > plotsetup(default); > plot([plot_yz, plot_yzLa, plot_yz1, plot_yz2], legend=["Real sol", > "Unwanted sol", "Backward Euler","Forward Euler"], colour=[green, > brown, red, black], style=[line,line,line, point]);

216


0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

0 –20 –40 –60 –80

Real sol Backward Euler Forward Euler

Now we can try h = 0.01. Here we find that forward Euler breaks down far quicker, inspite of having a step size.

217

APPENDIX B. MAPLE WORKSHEETS soln1h001x15 := GeneralBackwardEuler(u,f,ic,0.01,0.2): soln2h001x15 := GeneralForwardEuler(u,f,ic,0.01,0.2): plot_yz := []: plot_yzLa := []: plot_yz1 := []: plot_yz2 := []: for i from 1 by 1 to 21 do plot_yz := [op(plot_yz),[soln1h001x15[i,1], exp(-soln1h001x15 [i,1])]]: plot_yzLa := [op(plot_yzLa),[soln1h001x15[i,1], exp(-la*(soln1h001x15 [i,1]))]]: plot_yz1 := [op(plot_yz1),[soln1h001x15[i,1], soln1h001x15[i,2]]]: > plot_yz2 := [op(plot_yz2),[soln2h001x15[i,1],soln2h001x15[i,2]]]: > end do: > plotsetup("ps", plotoutput="h:\stiff-ploth001.eps", > plotoptions="color"); > plot([plot_yz, plot_yzLa, plot_yz1, plot_yz2], legend=["Real sol", > "Unwanted sol", "Backward Euler","Forward Euler"], colour=[green, > brown, red, black], style=[line,line,line, point]); > plotsetup(default); > plot([plot_yz, plot_yzLa, plot_yz1, plot_yz2], legend=["Real sol", > "Unwanted sol", "Backward Euler","Forward Euler"], colour=[green, > brown, red, black], style=[line,line,line, point]); > > > > > > > > > > > >

2

0.05

0.1

0.15

0 –2 –4 –6 –8


218

0.2


Now we can try h = .1e-2 Here we find that forward Euler is stable. This is not surprising as h < λ2 where 2 λ = .2e-2 > soln1h0001x15 := GeneralBackwardEuler(u,f,ic,0.001,10): > soln2h0001x15 := GeneralForwardEuler(u,f,ic,0.001,10): > plot_yz := []: > plot_yzLa := []: > plot_yz1 := []: > plot_yz2 := []: > for i from 1 by 100 to 10001 do > plot_yz := [op(plot_yz),[soln1h0001x15[i,1], exp(-soln1h0001x15 > [i,1])]]: > plot_yzLa := [op(plot_yzLa),[soln1h0001x15[i,1], > exp(-la*(soln1h0001x15 [i,1]))]]: > plot_yz1 := [op(plot_yz1),[soln1h0001x15[i,1], > soln1h0001x15[i,2]]]: > plot_yz2 := [op(plot_yz2),[soln2h0001x15[i,1],soln2h0001x15[i,2]]]: > end do: > #print(plot_yz); > plotsetup("ps", plotoutput="h:\stiff-ploth0001.eps", > plotoptions="color"); > plot([plot_yz, plot_yzLa, plot_yz1, plot_yz2], legend=["Real sol", > "Unwanted sol", "Backward Euler","Forward Euler"], colour=[green, > brown, red, black], style=[line,line,line, point]); > plotsetup(default); > plot([plot_yz, plot_yzLa, plot_yz1, plot_yz2], legend=["Real sol", > "Unwanted sol", "Backward Euler","Forward Euler"], colour=[green, > brown, red, black], style=[line,line,line, point]);

219


1 0.8 0.6 0.4 0.2 0

2

4

6

8

10


Next we can try λ = 106 We shall first let h = .1. > u := [x,y,z]; > f := [1,z,-la106*y - (la106 + 1)*z]; > ic := [0,1,-1]; u := [x, y, z] f := [1, z, −la106 y − (la106 + 1) z] ic := [0, 1, −1] >

la106 := 10^6; la106 := 1000000

If we limit the range for x = [0, .8], there are no sudden “jumps”.

220

APPENDIX B. MAPLE WORKSHEETS soln1Ah01x15 := GeneralBackwardEuler(u,f,ic,0.1,0.8): soln2Ah01x15 := GeneralForwardEuler(u,f,ic,0.1,0.8): plot_yz := []: plot_yzLa := []: plot_yz1 := []: plot_yz2 := []: for i from 1 by 1 to 9 do plot_yz := [op(plot_yz),[soln1Ah01x15[i,1], exp(-soln1Ah01x15 [i,1])]]: plot_yzLa := [op(plot_yzLa),[soln1Ah01x15[i,1], exp(-la*(soln1Ah01x15 [i,1]))]]: plot_yz1 := [op(plot_yz1),[soln1Ah01x15[i,1], soln1Ah01x15[i,2]]]: > plot_yz2 := [op(plot_yz2),[soln2Ah01x15[i,1],soln2Ah01x15[i,2]]]: > end do: > plot([plot_yz, plot_yzLa, plot_yz1, plot_yz2], legend=["Real sol", > "Unwanted sol", "Backward Euler","Forward Euler"], colour=[green, > brown, red, black], style=[line,line,line, point]); > > > > > > > > > > > >

1 0.9 0.8 0.7 0.6 0.5 0

0.2

0.4

0.6

0.8


However, when increasing the range to x = .9, we encounter the same problem seen earlier with forward Euler.

221

APPENDIX B. MAPLE WORKSHEETS soln1Ah01x16 := GeneralBackwardEuler(u,f,ic,0.1,0.9): soln2Ah01x16 := GeneralForwardEuler(u,f,ic,0.1,0.9): plot_yz := []: plot_yzLa := []: plot_yz1 := []: plot_yz2 := []: for i from 1 by 1 to 10 do plot_yz := [op(plot_yz),[soln1Ah01x16[i,1], exp(-soln1Ah01x16 [i,1])]]: plot_yzLa := [op(plot_yzLa),[soln1Ah01x16[i,1], exp(-la*(soln1Ah01x16 [i,1]))]]: plot_yz1 := [op(plot_yz1),[soln1Ah01x16[i,1], soln1Ah01x16[i,2]]]: > plot_yz2 := [op(plot_yz2),[soln2Ah01x16[i,1],soln2Ah01x16[i,2]]]: > end do: > plot([plot_yz, plot_yzLa, plot_yz1, plot_yz2], legend=["Real sol", > "Unwanted sol", "Backward Euler","Forward Euler"], colour=[green, > brown, red, black], style=[line,line,line, point]); > > > > > > > > > > > >

1000 800 600 400 200 0

0.2

0.4

0.6

0.8


Now we can try h = .1e-1 The solution generated by forward Euler very quickly becomes unstable.

222

APPENDIX B. MAPLE WORKSHEETS soln1Ah001x15 := GeneralBackwardEuler(u,f,ic,0.01,0.08): soln2Ah001x15 := GeneralForwardEuler(u,f,ic,0.01,0.08): plot_yz := []: plot_yzLa := []: plot_yz1 := []: plot_yz2 := []: for i from 1 by 1 to 9 do plot_yz := [op(plot_yz),[soln1Ah001x15[i,1], exp(-soln1Ah001x15 [i,1])]]: plot_yzLa := [op(plot_yzLa),[soln1Ah001x15[i,1], exp(-la*(soln1Ah001x15 [i,1]))]]: plot_yz1 := [op(plot_yz1),[soln1Ah001x15[i,1], soln1Ah001x15[i,2]]]: plot_yz2 := [op(plot_yz2),[soln2Ah001x15[i,1],soln2Ah001x15[i,2]]]: end do: plot([plot_yz, plot_yzLa, plot_yz1, plot_yz2], legend=["Real sol", > "Unwanted sol", "Backward Euler","Forward Euler"], colour=[green, > brown, red, black], style=[line,line,line, point]); > > > > > > > > > > > > > > > >

1 0.98 0.96 0.94 0

0.02

0.04

0.06

0.08


We could try a step size that satisfies the rule for stability say h = .1e-5 However, this is computationally very demanding, if we could run this, it would be stable. Note that h < λ2 where λ2 = .2e-5

223

APPENDIX B. MAPLE WORKSHEETS # procedure for running with a "stable" stepsize # in practice to takes too long to deliver a result stepsize := 0.000001: soln1Ah0001x15 := GeneralBackwardEuler(u,f,ic,stepsize,0.1): soln2Ah0001x15 := GeneralForwardEuler(u,f,ic,stepsize,0.1): plot_yz := []: plot_yzLa := []: plot_yz1 := []: plot_yz2 := []: for i from 1 by 100 to 10001 do plot_yz := [op(plot_yz),[soln1Ah0001x15[i,1], exp(-soln1Ah0001x15 [i,1])]]: plot_yzLa := [op(plot_yzLa),[soln1Ah0001x15[i,1], exp(-la*(soln1Ah0001x15 [i,1]))]]: plot_yz1 := [op(plot_yz1),[soln1Ah0001x15[i,1], soln1Ah0001x15[i,2]]]: plot_yz2 := [op(plot_yz2),[soln2Ah0001x15[i,1],soln2Ah0001x15[i,2]]]: end do: plot([plot_yz, plot_yzLa, plot_yz1, plot_yz2], legend=["Real sol", > "Unwanted sol", "Backward Euler","Forward Euler"], colour=[green, > brown, red, black], style=[line,line,line, point]); > > > > > > > > > > > > > > > > > > > >

B.9.4

How to Solve Analytically Using Maple

Comparing to the analytical solution Suppose we are trying to solve the ODE

(

d d2 y(x)) + (λ + 1) ( y(x)) + λ y(x) = 0 2 dx dx

d y(0) = −1. We specify the with the initial conditions given by y(0) = 1 and dx problem as follows: > ode := diff(diff(y(x), x), x) +(lambda + 1)*diff(y(x), x) + > lambda*y(x) = 0; 2

d d ode := ( dx 2 y(x)) + (λ + 1) ( dx y(x)) + λ y(x) = 0

We can use Maple to obtain an analytical solution for our ODE. This is achieved by using dsolve with different parameters: > ic1 := y(0) = 1; > ic2 := D(y)(0) = -1; > sol_true := dsolve({ode, ic1, ic2}); 224


ic1 := y(0) = 1 ic2 := D(y)(0) = −1 sol true := y(x) = e(−x) Appendix (required) Forward and Backward Euler methods > GeneralForwardEuler := proc(vars, f, init, h, interval) > local makeeqn, ans, solveEqns, assignList, diffValue, > newICSet,i,j,k, ic; > ic := init; > makeeqn := nops(u); > ans := []; > ans := [op(ans), ic]; > while op(1,ic) < interval do > solveEqns := []; > assignList := []; > for j from 1 by 1 to makeeqn do > assignList := [op(assignList), op(j, u) = op(j, ic)]; > end do; > for i from 1 by 1 to makeeqn do > diffValue := eval(op(i, f),assignList); > solveEqns := [op(solveEqns),op(i, u) = op(i, ic) + h * > (diffValue)]; > end do; > newICSet := solve(convert(solveEqns, set), convert(u, set)); > for k from 1 by 1 to makeeqn do > ic := subsop(k = eval(op(k, u), newICSet), ic); > end do; > ans := [op(ans), ic]; > end do; > return matrix(ans); > end proc:

225


GeneralBackwardEuler := proc(vars, f, init, h, interval) local makeeqn, ans, solveEqns, assignList, diffValue, newICSet,

i, > > > > > > > > > > > > > > > > > > > > > > > >

j, k, ic; ic := init; makeeqn := nops(u); ans := []; ans := [op(ans), ic]; while op(1,ic) < interval do solveEqns := []; assignList := []; for j from 1 by 1 to makeeqn do assignList := [op(assignList), op(j, u) = op(j, ic)]; end do; for i from 1 by 1 to makeeqn do diffValue := eval(op(i, f), assignList); solveEqns := [op(solveEqns), op(i, u) = op(i, ic) + h * (op(i, f))]; end do; newICSet := solve(convert(solveEqns, set),convert(u, set)); for k from 1 by 1 to makeeqn do ic := subsop(k = eval(op(k, u), newICSet), ic); end do; ans := [op(ans), ic]; end do; return matrix(ans); end proc:

226


B.10

6thorder.mws

6th Order Numerical Method S. Amen, P.A. Bilokon, A.J. Brinley Codd, M. Fofaria, T. Shah B.10.1

Initialization


B.10.2

Implementation

Implementation : 6th Order Method sixthOrder := proc(vars , f, init, h, interval) #Declaration of variables local nvars, ic, solveEqns, j,k,p, A,B,C,D,E,F,C_,D_,E_,F_,al, n_plus_half, n_minus_half, eq_plus_half,eq_minus_half,newICSet, variables,ans, eq1,eq2,eq3,eq4, eqn1,eqn2,eqn3,eqn4,eqn5, eqnA1,eqnA2,eqnA3,eqnA4,eqnA5,eqnA6, values, oldICSet; #What follows is the determination of the constants in the equations #provided (see report) eq1 := ((1/2+al)^2)/2 = C + D + E + F; 227


eq2 := ((1/2+al)^3)/6 = C + D*(1/2+al) + E*(1/2 - al); eq3 := ((1/2+al)^4)/24 = 1/2*(C + D*(1/2+al)^2 + E*(1/2 - al)^2); eq4 := ((1/2+al)^5)/120 = 1/6*(C + D*(1/2+al)^3 + E*(1/2 - al)^3); eqn1 := ((1/2-al)^2)/2 = C_ + D_ + E_ + F_; eqn2 := ((1/2-al)^3)/6 = C_ + D_*(1/2+al) + E_*(1/2 - al); eqn3 := ((1/2-al)^4)/24 = 1/2*(C_ + D_*(1/2+al)^2 + E_*(1/2 - al)^2); eqn4 := ((1/2-al)^5)/120 = 1/6*(C_ + D_*(1/2+al)^3 + E_*(1/2 - al)^3); eqnA1 := 1-2*A - 2*B = 0; eqnA2 := 1/2 - A - (B*(1/2+al)^1)/1 - (B*(1/2-al)^1)/1 = 0; eqnA3 := 1/6 - A/2 - (B*(1/2+al)^2)/2 - (B*(1/2-al)^2)/2 = 0; eqnA4 := 1/24 - A/6 - (B*(1/2+al)^3)/6 - (B*(1/2-al)^3)/6 = 0; eqnA5 := 1/120 - A/24 - (B*(1/2+al)^4)/24 - (B*(1/2-al)^4)/24 = 0; eqnA6 := 1/720 - A/120 - (B*(1/2+al)^5)/120 - (B*(1/2-al)^5)/120 = 0; values := evalf(solve( {eqnA1,eqnA2,eqnA3,eqnA4,eqnA5,eqnA6,eqn1,eqn2,eqn3,eqn4,eq1,eq2, eq3,eq4})); A := eval(A, values); B := eval(B, values); C := eval(C, values); D := eval(D, values); E := eval(E, values); F := eval(F, values); C_ := eval(C_, values); D_ := eval(D_, values); E_ := eval(E_, values); F_ := eval(F_, values); al := eval(al, values);

# determines the number of variables to solve for nvars := nops(vars); #some initialization ic := init; eq_plus_half := []; eq_minus_half := []; solveEqns := []; variables := []; ans := []; n_plus_half := [seq(plhlf[i],i=1..nvars)]; n_minus_half := [seq(minhlf[g],g=1..nvars)]; 228


oldICSet := {}; #more initialization, initialising the vectors of y_n+1/2+alpha and #y_n+1/2-alpha for k from 1 by 1 to nvars do eq_plus_half := [op(eq_plus_half),op(2,op(k,vars)) = op(k,n_plus_half)]; eq_minus_half := [op(eq_minus_half),op(2,op(k,vars)) = op(k,n_minus_half)]; oldICSet := oldICSet union {op(2,op(k,vars)) = op(2,op(k,ic))}; variables := [op(variables), op(2,op(k,vars))]; ans := [op(ans), [op(k,ic)]]; end do; #main iteration loop while op(1,op(1,ic)) < interval do #create simultaneous equations for y_n+1 etc for j from 1 by 1 to nvars do solveEqns := [op(solveEqns), op(2,op(j,vars)) = op(2,op(j,ic)) + h*op(3,op(j,ic)) + (h^2)*(B*(1/2-al)*subs(convert(eq_plus_half,set),op(3,op(j,f))) + B*(1/2+al)*subs(convert(eq_minus_half,set),op(3,op(j,f))) + A*eval(op(3,op(j,f)),oldICSet)), op(j,n_plus_half) = op(2,op(j,ic)) + (1/2 + al)*h*op(3,op(j,ic)) + (h^2)*(C*op(3,op(j,f))+D*subs(convert(eq_plus_half,set),op(3,op(j,f))) + E*subs(convert(eq_minus_half,set),op(3,op(j,f))) + F*eval(op(3,op(j,f)),oldICSet)), op(j,n_minus_half) = op(2,op(j,ic)) + (1/2 al)*h*op(3,op(j,ic)) + (h^2)*(C_*op(3,op(j,f))+ D_*subs(convert(eq_plus_half,set),op(3,op(j,f))) + E_*subs(convert(eq_minus_half,set),op(3,op(j,f))) + F_*eval(op(3,op(j,f)),oldICSet)) ]; 229


end do; #solve the equations we have created - NOTE fsolve is used because #it works a lot quicker than solve but the results are not different to the decimal #precision we used. newICSet := fsolve(convert(solveEqns,set),convert(n_plus_half,set) union convert(n_minus_half,set) union convert(variables,set));

#updating the variables to the new values for p from 1 by 1 to nvars do ic[p,3] := eval(op(3,op(p,ic)) + h*(A*op(3,op(p,f)) + B*subs(convert(eq_plus_half,set), op(3,op(p,f))) + B*subs(convert(eq_minus_half,set),op(3,op(p,f))) + A*eval(op(3,op(p,f)),oldICSet) ), newICSet ); ic[p,1] := op(1,op(p,ic)) + h; ic[p,2] := eval(op(2,op(p,vars)),newICSet); ans[p] := [op(ans[p]), op(p,ic)]; end do; #resetting values before next iteration solveEqns := []; oldICSet := newICSet; end do; return ans; end proc: Now to try the sixth order method on the following ODE : y’ = y , y ’ = -y, y(0) = 1, y (0) = 0 first to put it into a form admissible into our procedure:

> > > >

u := [[x,y,y_]]; f := [[1,y_,-1*y]]; ic := [[0,1,0]]; soln := sixthOrder(u,f,ic,0.1,10); 230


u := [[x, y, y ]] f := [[1, y , −y]] ic := [[0, 1, 0]] 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6

1 0.9950041656 0.9800665784 0.9553364899 0.9210609950 0.8775825632 0.8253356164 0.7648421888 0.6967067110 0.6216099700 0.5403023076 0.4535961231 0.3623577561 0.2674988303 0.1699671445 0.07073720312 -0.02919952105 -0.1288444933 -0.2272020940 -0.3232895665 -0.4161468366 -0.5048461050 -0.5885011181 -0.6662760226 -0.7373937172 -0.8011436175 -0.8568887558 -0.9040721447 -0.9422223436 -0.9709581685 -0.9899925003 -0.9991351535 -0.9982947792 -0.9874797734 -0.9667981962 -0.9364566909 -0.8967584200 231

0 -0.09983341664 -0.1986693308 -0.2955202067 -0.3894183424 -0.4794255388 -0.5646424737 -0.6442176877 -0.7173560915 -0.7833269104 -0.8414709857 -0.8912073611 -0.9320390872 -0.9635581868 -0.9854497316 -0.9974949884 -0.9995736050 -0.9916648125 -0.9738476330 -0.9463000898 -0.9092974290 -0.8632093688 -0.8084964058 -0.7457052140 -0.6754631822 -0.5984721457 -0.5155013732 -0.4273798814 -0.3349881510 -0.2392493298 -0.1411200083 -0.04158066230 0.05837414390 0.1577456950 0.2555411032 0.3507832292 0.4425204452


3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 6.0 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 7.0 7.1 7.2 7.3 7.4 7.5 7.6 7.7

-0.8481000354 -0.7909677156 -0.7259323077 -0.6536436243 -0.5748239498 -0.4902608243 -0.4007991747 -0.3073328723 -0.2107958014 -0.1121525284 -0.01238866442 0.08749898304 0.1865123696 0.2836621862 0.3779777440 0.4685166732 0.5543743387 0.6346928792 0.7086697781 0.7755658829 0.8347127898 0.8855195225 0.9274784368 0.9601702932 0.9832684453 0.9965421047 0.9998586445 0.9931849270 0.9765876340 0.9502326002 0.9143831564 0.8693974983 0.8157251077 0.7539022616 0.6845466732 0.6083513206 0.5260775228 0.4385473324 0.3466353219 0.2512598458 0.1533738644

232

0.5298361431 0.6118578935 0.6877661621 0.7568024986 0.8182771147 0.8715757764 0.9161659410 0.9516020784 0.9775301225 0.9936910087 0.9999232628 0.9961646144 0.9824526182 0.9589242802 0.9258146877 0.8834546610 0.8322674472 0.7727644922 0.7055403299 0.6312666418 0.5506855461 0.4646021824 0.3738766671 0.2794154999 0.1821625053 0.08308940319 -0.01681390089 -0.1165492061 -0.2151199902 -0.3115413665 -0.4048499244 -0.4941133557 -0.5784397697 -0.6569866047 -0.7289690467 -0.7936678710 -0.8504366284 -0.8987081041 -0.9379999855 -0.9679196811 -0.9881682431


> > > > > >

7.8 0.05395542200 -0.9985433547 7.9 -0.04600212515 -0.9989413513 8.0 -0.1455000343 -0.9893582561 8.1 -0.2435441553 -0.9698898203 8.2 -0.3391548635 -0.9407305660 8.3 -0.4313768485 -0.9021718428 8.4 -0.5192886586 -0.8545989167 8.5 -0.6020119083 -0.7984871208 8.6 -0.6787200539 -0.7343971055 8.7 -0.7486466531 -0.6629692371 8.8 -0.8110930225 -0.5849171991 8.9 -0.8654352184 -0.5010208618 9.0 -0.9111302718 -0.4121184896 9.1 -0.9477216127 -0.3190983657 9.2 -0.9748436326 -0.2228899164 9.3 -0.9922253370 -0.1244544247 9.4 -0.9996930540 -0.02477542552 9.5 -0.9971721683 0.07515112163 9.6 -0.9846878680 0.1743267836 9.7 -0.9623648920 0.2717606300 9.8 -0.9304262842 0.3664791340 9.9 -0.8891911645 0.4575358997 10.0 -0.8390715405 0.5440211180 #plot the phase plane plot then the error with the real solution plot_yz := []: err_plot_sixth := []: for i from 1 by 1 to 101 do plot_yz := [op(plot_yz),[soln[1,i,2],soln[1,i,3]]]: err_plot_sixth := [op(err_plot_sixth),[soln[1,i,1],abs(soln[1,i,2]

> > > >

cos(soln[1,i,1]))]]: end do: plot(plot_yz, labels=["y", "y’"]); plot(err_plot_sixth,labels=["x", "error"]);

233


1

y’

–1

0.5

–0.5

0.5 y

1

–0.5

–1

1.2e–08 1e–08 8e–09 error 6e–09 4e–09 2e–09 0

2

4

6 x

234

8

10

APPENDIX B. MAPLE WORKSHEETS #Now solve Kepler’s equations for movement of sun around Earth #Note we use a stepsize of 0.1 here but we have also calculated with a > #stepsize of 0.01 but do not include it here due to space constraints. > u := [[x,y,y_],[x,z,z_]]; > f := [[1,y_,-1*4*Pi^2*y/(y^2 + z^2)^(3/2)],[1,z_,-1*4*Pi^2*z/(z^2 + > y^2)^(3/2)]]; > ic := [[0,0,2*Pi],[0,1,-1*Pi]]; > soln2 := sixthOrder(u,f,ic,0.1,10); u := [[x, y, y ], [x, z, z ]] > >

f := [[1, y , −

4 π2 y 4 π2 z ], [1, z , − ]] (y 2 + z 2 )(3/2) (y 2 + z 2 )(3/2)

ic := [[0, 0, 2 π], [0, 1, −π]] soln2 := [[0, 0, 2*Pi], [.1, .5605799939, 2*Pi-.2374104392*Pi^2], [.2, .5644782548, 2*Pi-1.029627099*Pi^2], [.3, .6653334426e-2, 2*Pi-1.272781521*Pi^2], [.4, -.5949377183, 2*Pi-1.202015944*Pi^2], [.5, -1.096332140, 2*Pi-1.085886587*Pi^2], [.6, -1.482294470, 2*Pi-.9704724664*Pi^2], [.7, -1.757514234, 2*Pi-.8615296079*Pi^2], [.8, -1.927999161, 2*Pi-.7578808130*Pi^2], [.9, -1.997844072, 2*Pi-.6571792502*Pi^2], [1.0, -1.968738125, 2*Pi-.5569525968*Pi^2], [1.1, -1.839929424, 2*Pi-.4547255302*Pi^2], [1.2, -1.608225155, 2*Pi-.3480464866*Pi^2], [1.3, -1.268178990, 2*Pi-.2349751997*Pi^2], [1.4, -.8136410995, 2*Pi-.1169827847*Pi^2], [1.5, -.2468823068, 2*Pi-.1514189210e-1*Pi^2], [1.6, .3674082279, 2*Pi-.6775558585e-1*Pi^2], [1.7, .6624812407, 2*Pi-.7259298174*Pi^2], [1.8, .2513570909, 2*Pi-1.246717527*Pi^2], [1.9, -.3672584132, 2*Pi-1.242135037*Pi^2], [2.0, -.9131431529, 2*Pi-1.132695990*Pi^2], [2.1, -1.344369231, 2*Pi-1.014980786*Pi^2], [2.2, -1.662020521, 2*Pi-.9030881391*Pi^2], [2.3, -1.872379246, 2*Pi-.7972865015*Pi^2], [2.4, -1.980466214, 2*Pi-.6954369625*Pi^2], [2.5, -1.988966499, 2*Pi-.5950676657*Pi^2], [2.6, -1.898117058, 2*Pi-.4936992155*Pi^2], 235


[2.7, [2.8, [2.9, [3.0, [3.1, [3.2, [3.3, [3.4, [3.5, [3.6, [3.7, [3.8, [3.9, [4.0, [4.1, [4.2, [4.3, [4.4, [4.5, [4.6, [4.7, [4.8, [4.9, [5.0, [5.1, [5.2, [5.3, [5.4, [5.5, [5.6, [5.7, [5.8, [5.9, [6.0, [6.1, [6.2, [6.3, [6.4, [6.5, [6.6, [6.7, [6.8,

-1.705708978, 2*Pi-.3888584200*Pi^2], -1.407160586, 2*Pi-.2782982392*Pi^2], -.9962903398, 2*Pi-.1615342521*Pi^2], -.4700742039, 2*Pi-.4824854347e-1*Pi^2], .1438324164, 2*Pi-.8037857570e-2*Pi^2], .6336353554, 2*Pi-.3981561885*Pi^2], .4613555201, 2*Pi-1.146343012*Pi^2], -.1380933150, 2*Pi-1.268416572*Pi^2], -.7217771751, 2*Pi-1.176749008*Pi^2], -1.196909161, 2*Pi-1.058989368*Pi^2], -1.556899677, 2*Pi-.9448606717*Pi^2], -1.807550640, 2*Pi-.8373039445*Pi^2], -1.954645252, 2*Pi-.7346278581*Pi^2], -2.001819294, 2*Pi-.6344092359*Pi^2], -1.950276618, 2*Pi-.5341749008*Pi^2], -1.798783774, 2*Pi-.4314562525*Pi^2], -1.543687496, 2*Pi-.3238537984*Pi^2], -1.179238202, 2*Pi-.2097727572*Pi^2], -.6998839932, 2*Pi-.9257925803e-1*Pi^2], -.1135170122, 2*Pi-.4203319996e-2*Pi^2], .4800749223, 2*Pi-.1434952173*Pi^2], .6242670188, 2*Pi-.9044337678*Pi^2], .1198407323, 2*Pi-1.268578366*Pi^2], -.4932055554, 2*Pi-1.221723399*Pi^2], -1.015489269, 2*Pi-1.107332083*Pi^2], -1.421959967, 2*Pi-.9904610301*Pi^2], -1.716118727, 2*Pi-.8799750056*Pi^2], -1.904221392, 2*Pi-.7752186676*Pi^2], -1.990835803, 2*Pi-.6739056253*Pi^2], -1.978136939, 2*Pi-.5735544172*Pi^2], -1.865851072, 2*Pi-.4716869838*Pi^2], -1.651262327, 2*Pi-.3658440094*Pi^2], -1.329346628, 2*Pi-.2539387822*Pi^2], -.8939341549, 2*Pi-.1364642042*Pi^2], -.3445907445, 2*Pi-.2862277887e-1*Pi^2], .2731938609, 2*Pi-.3316682164e-1*Pi^2], .6650751167, 2*Pi-.5790086680*Pi^2], .3469040908, 2*Pi-1.215638587*Pi^2], -.2699593345, 2*Pi-1.256292476*Pi^2], -.8335746530, 2*Pi-1.152285248*Pi^2], -1.283952584, 2*Pi-1.033916358*Pi^2], -1.619678151, 2*Pi-.9208058036*Pi^2], 236


[6.9, -1.847057395, 2*Pi-.8141370002*Pi^2], [7.0, -1.971545179, 2*Pi-.7119030540*Pi^2], [7.1, -1.996307329, 2*Pi-.6116356749*Pi^2], [7.2, -1.922060097, 2*Pi-.5108563388*Pi^2], [7.3, -1.747079675, 2*Pi-.4070987113*Pi^2], [7.4, -1.467252725, 2*Pi-.2980483262*Pi^2], [7.5, -1.076625669, 2*Pi-.1826683812*Pi^2], [7.6, -.5708883379, 2*Pi-.6736665542e-1*Pi^2], [7.7, .3268095960e-1, 2*Pi-.1754572195e-2*Pi^2], [7.8, .5795569555, 2*Pi-.2690989753*Pi^2], [7.9, .5442662449, 2*Pi-1.059680228*Pi^2], [8.0, -.2615007433e-1, 2*Pi-1.274012817*Pi^2], [8.1, -.6251321810, 2*Pi-1.197173344*Pi^2], [8.2, -1.121081457, 2*Pi-1.080001613*Pi^2], [8.3, -1.501161124, 2*Pi-.9644985199*Pi^2], [8.4, -1.770527927, 2*Pi-.8556529475*Pi^2], [8.5, -1.935269812, 2*Pi-.7521171522*Pi^2], [8.6, -1.999472103, 2*Pi-.6515031892*Pi^2], [8.7, -1.964793651, 2*Pi-.5513309684*Pi^2], [8.8, -1.830451521, 2*Pi-.4491306564*Pi^2], [8.9, -1.593235545, 2*Pi-.3424762620*Pi^2], [9.0, -1.247734027, 2*Pi-.2295226079*Pi^2], [9.1, -.7880366725, 2*Pi-.1121344125*Pi^2], [9.2, -.2175973750, 2*Pi-.1333049459e-1*Pi^2], [9.3, .3926241938, 2*Pi-.8188096132e-1*Pi^2], [9.4, .6576682978, 2*Pi-.7651764285*Pi^2], [9.5, .2234966729, 2*Pi-1.254351920*Pi^2], [9.6, -.3958792047, 2*Pi-1.238957304*Pi^2], [9.7, -.9373985160, 2*Pi-1.127558567*Pi^2], [9.8, -1.363339047, 2*Pi-1.009511312*Pi^2], [9.9, -1.675577397, 2*Pi-.8976084412*Pi^2], [10.0, -1.880546995, 2*Pi-.7918456882*Pi^2]], [[0, [.1, [.2, [.3, [.4, [.5, [.6, [.7, [.8,

1, -Pi], .4512244308, -Pi-.4959786788*Pi^2], -.4437144412, -Pi-.5005407941*Pi^2], -.9973464448, -Pi-.422715980e-2*Pi^2], -1.153922018, -Pi+.2917477049*Pi^2], -1.094154972, -Pi+.4506172138*Pi^2], -.9146167716, -Pi+.5417965370*Pi^2], -.6652030097, -Pi+.5954112560*Pi^2], -.3755066960, -Pi+.6248898790*Pi^2], 237


[.9, -.6594681436e-1, -Pi+.6362848867*Pi^2], [1.0, .2469185786, -Pi+.6316826950*Pi^2], [1.1, .5471158182, -Pi+.6102249333*Pi^2], [1.2, .8162604593, -Pi+.5676961972*Pi^2], [1.3, 1.029347741, -Pi+.4943011044*Pi^2], [1.4, 1.146668564, -Pi+.3684165993*Pi^2], [1.5, 1.095194277, -Pi+.1400105723*Pi^2], [1.6, .7283543005, -Pi-.2866884376*Pi^2], [1.7, -.9386788847e-1, -Pi-.6299477800*Pi^2], [1.8, -.8370609858, -Pi-.1826729322*Pi^2], [1.9, -1.124284657, -Pi+.1980818695*Pi^2], [2.0, -1.133219377, -Pi+.3998455252*Pi^2], [2.1, -.9920256065, -Pi+.5126561156*Pi^2], [2.2, -.7647316164, -Pi+.5787399571*Pi^2], [2.3, -.4871715589, -Pi+.6165105033*Pi^2], [2.4, -.1826189828, -Pi+.6343343494*Pi^2], [2.5, .1312352199, -Pi+.6356424344*Pi^2], [2.6, .4384526233, -Pi+.6206900854*Pi^2], [2.7, .7218857584, -Pi+.5866801188*Pi^2], [2.8, .9596264929, -Pi+.5263615274*Pi^2], [2.9, 1.118855396, -Pi+.4237659910*Pi^2], [3.0, 1.142452045, -Pi+.2426429051*Pi^2], [3.1, .9174780817, -Pi-.9818486116e-1*Pi^2], [3.2, .2565888697, -Pi-.5897939440*Pi^2], [3.3, -.6146985965, -Pi-.3824276333*Pi^2], [3.4, -1.059285967, -Pi+.8201129865e-1*Pi^2], [3.5, -1.152996255, -Pi+.3375063628*Pi^2], [3.6, -1.058989831, -Pi+.4764985709*Pi^2], [3.7, -.8596470118, -Pi+.5570177709*Pi^2], [3.8, -.5986823060, -Pi+.6040427023*Pi^2], [3.9, -.3028886775, -Pi+.6288199687*Pi^2], [4.0, .8593783536e-2, -Pi+.6363224077*Pi^2], [4.1, .3196491828, -Pi+.6279460542*Pi^2], [4.2, .6141091639, -Pi+.6021849174*Pi^2], [4.3, .8726529092, -Pi+.5539057957*Pi^2], [4.4, 1.068087772, -Pi+.4715550191*Pi^2], [4.5, 1.155834682, -Pi+.3294569724*Pi^2], [4.6, 1.051572760, -Pi+.6804147613e-1*Pi^2], [4.7, .5894164639, -Pi-.4023250794*Pi^2], [4.8, -.2888880795, -Pi-.5775045772*Pi^2], [4.9, -.9316775324, -Pi-.8084008120e-1*Pi^2], [5.0, -1.143519583, -Pi+.2524949572*Pi^2], 238


[5.1, [5.2, [5.3, [5.4, [5.5, [5.6, [5.7, [5.8, [5.9, [6.0, [6.1, [6.2, [6.3, [6.4, [6.5, [6.6, [6.7, [6.8, [6.9, [7.0, [7.1, [7.2, [7.3, [7.4, [7.5, [7.6, [7.7, [7.8, [7.9, [8.0, [8.1, [8.2, [8.3, [8.4, [8.5, [8.6, [8.7, [8.8, [8.9, [9.0, [9.1, [9.2,

-1.112410709, -Pi+.4295782443*Pi^2], -.9486353343, -Pi+.5299550342*Pi^2], -.7080879494, -Pi+.5888624613*Pi^2], -.4230318160, -Pi+.6218384093*Pi^2], -.1151152613, -Pi+.6359277819*Pi^2], .1986185220, -Pi+.6338043744*Pi^2], .5022272587, -Pi+.6151093010*Pi^2], .7779051565, -Pi+.5762818504*Pi^2], 1.002096411, -Pi+.5087297465*Pi^2], 1.138354222, -Pi+.3935537158*Pi^2], 1.121277299, -Pi+.1873831861*Pi^2], .8199168239, -Pi-.2008580031*Pi^2], .6142804843e-1, -Pi-.6335323644*Pi^2], -.7484432935, -Pi-.2675956188*Pi^2], -1.099863951, -Pi+.1518542295*Pi^2], -1.142336746, -Pi+.3753605379*Pi^2], -1.019563206, -Pi+.4986519321*Pi^2], -.8030149510, -Pi+.5704673819*Pi^2], -.5317133484, -Pi+.6118748512*Pi^2], -.2303494078, -Pi+.6324179100*Pi^2], .8282210742e-1, -Pi+.6361719580*Pi^2], .3918398474, -Pi+.6238885178*Pi^2], .6800492239, -Pi+.5933597342*Pi^2], .9267849213, -Pi+.5383379024*Pi^2], 1.101931618, -Pi+.4449985061*Pi^2], 1.154461074, -Pi+.2822940785*Pi^2], .9852957363, -Pi-.2099682400e-1*Pi^2], .4120886512, -Pi-.5188006246*Pi^2], -.4819342769, -Pi-.4767448263*Pi^2], -1.011175171, -Pi+.1637749707e-1*Pi^2], -1.152209765, -Pi+.3035061864*Pi^2], -1.083381691, -Pi+.4577014550*Pi^2], -.8982814277, -Pi+.5461959422*Pi^2], -.6454556873, -Pi+.5980258184*Pi^2], -.3538717957, -Pi+.6261480024*Pi^2], -.4365344015e-1, -Pi+.6363795234*Pi^2], .2687546077, -Pi+.6306578619*Pi^2], .5673572920, -Pi+.6079912873*Pi^2], .8336013331, -Pi+.5639877188*Pi^2], 1.042080322, -Pi+.4885329677*Pi^2], 1.152210297, -Pi+.3593024227*Pi^2], 1.088995996, -Pi+.1246554165*Pi^2], 239


[9.3, .7024176002, -Pi-.3106833270*Pi^2], [9.4, -.1340893897, -Pi-.6234253252*Pi^2], [9.5, -.8574700739, -Pi-.1601212513*Pi^2], [9.6, -1.127182380, -Pi+.2113890134*Pi^2], [9.7, -1.125934376, -Pi+.4077614331*Pi^2], [9.8, -.9785084571, -Pi+.5176286832*Pi^2], [9.9, -.7472831966, -Pi+.5818513852*Pi^2], [10.0, -.4673512590, -Pi+.6182606986*Pi^2]]] > #Plotting y against z > #notice the graph here is not a clear curve. But when we used > #a stepsize of 0.01 we got a well defined curve (see sixth order > #method chapter in report). > plot_yz := []: > for i from 1 by 1 to 101 do > plot_yz := [op(plot_yz),[soln2[2,i,2],soln2[1,i,2]]]: > end do: > plot(plot_yz,labels=["Z", "Y"]);

0.5 –1

Z 0.5

–0.5 0 –0.5 Y

–1 –1.5 –2

240

1

APPENDIX B. MAPLE WORKSHEETS #This is the solution to the Kepler’s Equation problem posed by professor Cash #Note we use a stepsize of 0.1 here but we have also calculated with a > #stepsize of 0.01 but do not include it here due to space constraints. > u := [[x,y,y_],[x,z,z_]]; > f := [[1,y_,-1*y/(y^2 + z^2)^(3/2)],[1,z_,-1*z/(z^2 + y^2)^(3/2)]]; > ic := [[0,1-0.6,0],[0,0,((1+0.6)/(1-0.6))^(1/2)]]; > soln3 := sixthOrder(u,f,ic,0.1,10); u := [[x, y, y ], [x, z, z ]] y z f := [[1, y , − 2 ], [1, z , − 2 ]] (3/2) 2 (y + z ) (y + z 2 )(3/2) > > >

ic := [[0, 0.4, 0], [0, 0, 2.000000000]] soln3 := [[[0, .4, 0], [.1, .3698264891, -.5830974346], [.2, .2898770734, -.9788160158], [.3, .1807788640, -1.175483264], [.4, .5902691976e-1, -1.244028812], [.5, -.6577535598e-1, -1.244129868], [.6, -.1886807494, -1.210145482], [.7, -.3072723729, -1.159920919], [.8, -.4204308881, -1.102534655], [.9, -.5276972046, -1.042613764], [1.0, -.6289489894, -.9825160250], [1.1, -.7242345070, -.9234214177], [1.2, -.8136856375, -.8658871153], [1.3, -.8974715796, -.8101365757], [1.4, -.9757739797, -.7562138084], [1.5, -1.048773606, -.7040674264], [1.6, -1.116643360, -.6535971580], [1.7, -1.179544811, -.6046798039], [1.8, -1.237626696, -.5571837117], [1.9, -1.291024518, -.5109767292], [2.0, -1.339860724, -.4659304122], [2.1, -1.384245192, -.4219220663], [2.2, -1.424275857, -.3788355378], [2.3, -1.460039369, -.3365612859], [2.4, -1.491611742, -.2949960533], [2.5, -1.519058957, -.2540423189], [2.6, -1.542437503, -.2136076420], 241


[2.7, [2.8, [2.9, [3.0, [3.1, [3.2, [3.3, [3.4, [3.5, [3.6, [3.7, [3.8, [3.9, [4.0, [4.1, [4.2, [4.3, [4.4, [4.5, [4.6, [4.7, [4.8, [4.9, [5.0, [5.1, [5.2, [5.3, [5.4, [5.5, [5.6, [5.7, [5.8, [5.9, [6.0, [6.1, [6.2, [6.3, [6.4, [6.5, [6.6, [6.7, [6.8,

-1.561794848, -.1736039606], -1.577169845, -.1339468792], -1.588593061, -.9455496221e-1], -1.596087036, -.5534904207e-1], -1.599666482, -.1625154259e-1], -1.599338403, .2281418488e-1], -1.595102159, .6192451745e-1], -1.586949456, .1011562135], -1.574864280, .1405870698], -1.558822756, .1802965869], -1.538792947, .2203666486], -1.514734580, .2608822165], -1.486598700, .3019320379], -1.454327261, .3436093582], -1.417852637, .3860126177], -1.377097070, .4292460947], -1.331972063, .4734204249], -1.282377722, .5186528800], -1.228202099, .5650672002], -1.169320583, .6127926393], -1.105595452, .6619616361], -1.036875763, .7127051121], -.9629978997, .7651436607], -.8837873272, .8193715716], -.7990625149, .8754282092], -.7086427458, .9332466827], -.6123629276, .9925608856], -.5101011940, 1.052734338], -.4018303277, 1.112438191], -.2877146044, 1.169030477], -.1682953837, 1.217330862], -.4485346768e-1, 1.247168884], .7987700894e-1, 1.238571398], .2003406708, 1.153424813], .3058956236, .9277881792], .3788897678, .4952863805], .3991194705, -.1047499349], .3593366706, -.6650924545], .2730504559, -1.024068658], .1608769623, -1.193980660], .3810113392e-1, -1.247617556], -.8663831348e-1, -1.240166002], 242


[6.9, -.2089417001, -1.202520883], [7.0, -.3266755180, -1.150618844], [7.1, -.4388645638, -1.092564441], [7.2, -.5451232317, -1.032490685], [7.3, -.6453666116, -.9725007552], [7.4, -.7396616437, -.9136408398], [7.5, -.8281489574, -.8563966810], [7.6, -.9110011776, -.8009533014], [7.7, -.9884005154, -.7473341296], [7.8, -1.060526784, -.6954770315], [7.9, -1.127551159, -.6452764660], [8.0, -1.189633122, -.5966070253], [8.1, -1.246919196, -.5493365447], [8.2, -1.299542666, -.5033332738], [8.3, -1.347623830, -.4584696332], [8.4, -1.391270517, -.4146239959], [8.5, -1.430578720, -.3716813272], [8.6, -1.465633259, -.3295331721], [8.7, -1.496508423, -.2880772784], [8.8, -1.523268564, -.2472170254], [8.9, -1.545968624, -.2068607573], [9.0, -1.564654598, -.1669210776], [9.1, -1.579363927, -.1273141365], [9.2, -1.590125814, -.8795892715e-1], [9.3, -1.596961474, -.4877659618e-1], [9.4, -1.599884316, -.968976918e-2], [9.5, -1.598900059, .2937811134e-1], [9.6, -1.594006779, .6850344074e-1], [9.7, -1.585194894, .1077631059], [9.8, -1.572447081, .1472351426], [9.9, -1.555738128, .1869994038], [10.0, -1.535034719, .2271382411]], [[0, [.1, [.2, [.3, [.4, [.5, [.6, [.7, [.8,

0, 2.000000000], .1950361828, 1.855667313], .3649597191, 1.527446218], .4998455628, 1.175137257], .6016951835, .8720416440], .6762738355, .6289954578], .7291929672, .4368728508], .7649566910, .2840775736], .7869964056, .1610032350], 243


[.9, .7979064339, .6046678468e-1], [1.0, .7996651134, -.2276224050e-1], [1.1, .7938028752, -.9249145101e-1], [1.2, .7815227326, -.1515199564], [1.3, .7637854944, -.2019334272], [1.4, .7413702573, -.2453085000], [1.5, .7149178202, -.2828539216], [1.6, .6849622743, -.3155086880], [1.7, .6519543290, -.3440110608], [1.8, .6162787908, -.3689477404], [1.9, .5782678499, -.3907893734], [2.0, .5382113207, -.4099165269], [2.1, .4963646408, -.4266389259], [2.2, .4529552009, -.4412098588], [2.3, .4081874186, -.4538370735], [2.4, .3622468587, -.4646910879], [2.5, .3153036233, -.4739115685], [2.6, .2675151804, -.4816122458], [2.7, .2190287600, -.4878847036], [2.8, .1699834166, -.4928012864], [2.9, .1205118378, -.4964173041], [3.0, .7074196202e-1, -.4987726601], [3.1, .2079845813e-1, -.4998929954], [3.2, -.2919588766e-1, -.4997904098], [3.3, -.7911883762e-1, -.4984637982], [3.4, -.1288473835, -.4958988204], [3.5, -.1782564063, -.4920675012], [3.6, -.2272172764, -.4869274421], [3.7, -.2755963545, -.4804206004], [3.8, -.3232533454, -.4724715715], [3.9, -.3700394518, -.4629852768], [4.0, -.4157952616, -.4518439219], [4.1, -.4603482871, -.4389030374], [4.2, -.5035100513, -.4239863445], [4.3, -.5450725893, -.4068790887], [4.4, -.5848041842, -.3873193469], [4.5, -.6224441069, -.3649866158], [4.6, -.6576960379, -.3394866997], [4.7, -.6902197374, -.3103314955], [4.8, -.7196203566, -.2769116418], [4.9, -.7454345393, -.2384590503], [5.0, -.7671120929, -.1939949017], 244


[5.1, [5.2, [5.3, [5.4, [5.5, [5.6, [5.7, [5.8, [5.9, [6.0, [6.1, [6.2, [6.3, [6.4, [6.5, [6.6, [6.7, [6.8, [6.9, [7.0, [7.1, [7.2, [7.3, [7.4, [7.5, [7.6, [7.7, [7.8, [7.9, [8.0, [8.1, [8.2, [8.3, [8.4, [8.5, [8.6, [8.7, [8.8, [8.9, [9.0, [9.1, [9.2,

-.7839914696, -.1422564900], -.7952664771, -.8159298257e-1], -.7999404197, -.981526657e-2], -.7967620635, .7602159441e-1], -.7841353162, .1799316545], -.7599915991, .3074340444], -.7216127723, .4660937243], -.6654026130, .6659320022], -.5866613066, .9186392812], -.4796365805, 1.231778514], -.3388007210, 1.587681249], -.1635108539, 1.897689980], .3356385267e-1, 1.995603434], .2258101951, 1.808374776], .3901027883, 1.466791250], .5191159120, 1.120021842], .6159610339, .8271742013], .6865368678, .5934973222], .7362917267, .4087556551], .7695365385, .2615536243], .7895449088, .1427060209], .7987936345, .4539690902e-1], .7991755309, -.3533072483e-1], .7921587362, -.1030897102], .7789004221, -.1605415789], .7603271993, -.2096742593], .7371922816, -.2519941275], .7101166020, -.2886587077], .6796188028, -.3205690980], .6461374358, -.3484351068], .6100476374, -.3728217655], .5716738352, -.3941829903], .5312995628, -.4128862585], .4891751438, -.4292309197], .4455237852, -.4434619329], .4005464731, -.4557802725], .3544259555, -.4663508741], .3073300282, -.4753087391], .2594142826, -.4827636406], .2108244390, -.4888037490], .1616983624, -.4934984110], .1121678341, -.4969002494], 245


[9.3, .6236014301e-1, -.4990467044], [9.4, .1239954686e-1, -.4999611030], [9.5, -.3759135260e-1, -.4996533138], [9.6, -.8749028049e-1, -.4981200212], [9.7, -.1371739689, -.4953446329], [9.8, -.1865168091, -.4912968177], [9.9, -.2353894382, -.4859316474], [10.0, -.2836572188, -.4791882982]]] > plot_yz := []: > for i from 1 by 1 to 101 do > plot_yz := [op(plot_yz),[soln3[2,i,2],soln3[1,i,2]]]: > end do: > plot(plot_yz,labels=["Z", "Y"] );

–0.8 –0.6 –0.4 –0.2

0.2 0

–0.5 Y –1

–1.5

246

Z 0.4

0.6

0.8

Bibliography [1] H. G. Press B. P. Flannery S. A. Teukolsky W. T. Vetterling. Numerical Recipes in C: The Art of Scientific Computing. Cambridge University Press, 2nd edition, 1992. [2] E. J. Kostelich D. Armbruster. Introductory Differential Equations: From Linearity to Chaos. Addison Wesley Longman, 1996. [3] K. A. Stroud D. J. Booth. Engineering Mathematics. PALGRAVE, 5th edition, 2001. [4] D. W. Jordan P. Smith. Nonlinear Ordinary Differential Equations. Clarendon Press, Oxford, 2nd edition, 1987. [5] H. G. Press S. A. Teukolsky S. A. Vetterling B. P. Flannery. Numerical Recipes in C++. Cambridge University Press, 2nd edition, 2002. [6] H.G. Press S. A. Teukolsky W. T. Vetterling B. P. Flannery. Numerical Recipes in Fortran 90. Cambridge University Press, 2nd edition, 1996. [7] R. Bronson. Schaum’s Outline of Theory and Problems of Differential Equations. Schaum’s Outlines. McGraw-Hill, 2nd edition, 1993. [8] J. R. Cash. A variable step runge-kutta-nystr¨ om integrator for reversible systems of second order initial value problems. 2003. [9] E. Kreyszig. Advanced Engineering Mathematics. Wiley Text Books, 8th edition, 1998. [10] J. D. Lambert. Computational Methods in Ordinary Differential Equations. John Wiley Sons, 1973. [11] R. Redheffer. Differential Equations: Theory and Applications. Jones Bartlett Publishers, 1991. [12] G. Stephenson. Mathematical Methods for Science Students. Longman, 2nd edition, 1973.

247

Index implementation of numerical method sixth order, 47 implicit numerical method, 14

absolute error, 15 Backward Euler’s method, 13 circle problem, 47

Kepler’s Equations, 6, 31, 32, 36, 49, 50 Kepler, Johannes, 55 Kutta, Martin Wilhelm, 55

degree of differential equation, 7 differential equation degree of, 7 higher order, 7 linear, 7 non-linear, 7 order, 7 ordinary, 7 partial, 7 system, 7 differential equations first order, 27, 31, 39 second order system, 27

local error, 15 Newton, Sir Isaac, 51 numerical method, 11 asymmetric, 19 Backward Euler’s, 13 BackwardEuler, 40 Euler’s, see Frward Euler’s12 explicit, 14 Forward Euler, 40 Forward Euler’s, 12 implicit, 14 order of, 15 symmetric, 18 Trapezium Rule, 15 Trapezoidal Rule, see Trapezium Rule

error absolute, 15 global, 15 local, 15, 43, 44 roundoff, 15 Euler’s method, see Frward Euler’s method12 Euler, Leonhard, 52 explicit numerical method, 14

order of differential equation, 7 of numerical method, 15, 43 sixth order, 43

Forward Euler, 40 Forward Euler’s method, 12

phase plane plot, 8

global error, 15

roundoff error, 15 Runge, Carle David Tolmé, 54 248

INDEX

INDEX

Simpson’s Rule, 28 stiff equation, 39 symmetric, 28 symmetricity, 43 symmetry, 18 system of differential equations, 7 Taylor series, 44–46 time, 8 time series plot, 8 Trapeizum Rule, 35 Trapezium Rule, 15 Trapezoidal Rule, see Trapezium Rule

249