Els apunts del campus virtual, Apuntes de Ingeniería de Sistemas Audiovisuales

¡Descarga Els apunts del campus virtual y más Apuntes en PDF de Ingeniería de Sistemas Audiovisuales solo en Docsity!

Numerical Methods

Josep Blat

Universitat Pompeu Fabra

Contents

Chapter 1

Introduction and Basic

Concepts

1.1 Basic Concepts

Numerical Methods (NM) is a part of Mathematics mostly devoted to obtain solutions of equations in a practical way. We recall that an equation usually involves one or more unknowns, and that the value(s) of the unknown(s) that satisfy the equation are called solutions. One of the concerns of NM is to obtain solutions efficiently, i.e., by using algorithms (procedures to compute the solution) which are as fast as possible. Another important concern is to estimate the errors, i.e., how much the numerical solution differs from the exact solution of the equation. A source of errors which always exists is due to the fact that most real numbers can only be represented in an approximate way in the computer. This type of errors are called rounding errors. Other errors can come from the computer operations themselves, which are sometimes based on approximate algorithms. The algorithms can be another source of errors. These errors derived from the algorithms are called discretization errors. We shall study the so called convergence and stability of the algorithms. A final important aspect is the quality of the computer code. Fundamental properties are reliability and robustness; good programs should be portable and maintainable as well.

1.2 Orientation of the Course

The focus of this course is on numerical methods for Differential Equations (DEs), as their applications are an important motivation. The most interesting DEs stem from modelling physical phenomena. An important current applica- tion of this type of DEs is visual simulation. They are frequently used to generate the visual effects appearing in films and games. Some examples that we will consider in this course are the simulation of water and smoke. The vi- sual simulation should also help to understand better the errors, as well as the convergence and stability concepts.

CHAPTER 1. INTRODUCTION AND BASIC CONCEPTS 4

The consideration of efficiency that we mentioned is especially important, for instance, in games, where the simulation has to be real-time, i.e., it has to produce around 30 frames per second, which means that a result has to be obtained in 1/30th of a second. In order to understand better the applications, we shall discuss the physical models and the derivation of the corresponding DEs, and, if convenient, their exact solutions, also called analytical solutions, which we distinguish from their approximate or numerical solutions. If the DEs include partial derivatives they are called Partial Differential Equations (PDEs); if not, they are called Ordinary Differential Equations (ODEs). The following image represents the different fields involved in this process.

1.3 Basic ODE Concepts

Definition 1. Let x be a function of one variable, usually denoted as t, which is n times

differentiable. An expression such as F (t, x, x′,... , x(n)) = 0, n ≥ 1 , is calle an (ordinary) differential equation and n is the order of the equation. 2

A first example of first order ODE obtained from a physical model that rep- resenting the uniform movement: x′(t) = v, if v is constant. The ODE repre- senting the uniformly accelerated movement is second order: x′′(t) = a, a being constant. A function f (t) such that F (t, f (t), f ′(t),... , f (n)(t)) = 0 for every t is a solution of the equation. Remark that the solution of a DE is a function, when the solution of an algebraic equation is usually just a number (or several of them, one for each unknown). We call this solution exact and sometimes analytical. The numerical/approximate solution calculated with a computer through a numerical method cannot be obtained for all the infinite values of t: we can only obtain the solution for a finite number of values of t, usually tn = n∆t; and the solution is only an approximation of the exact solution. In the expressions used, we are implicitly assuming that a solution of the ODE exists and, furthermore, that it is unique. Indeed, this is the case for ODEs coming from physical models and under suitable conditions. Mathematically, Picard’s theorem shows that the solution of a first order ODE with a given initial value x 0 = x(0) exists and is unique, for t ∈ (0, a). A successive

CHAPTER 1. INTRODUCTION AND BASIC CONCEPTS 6

fact unique, as the particle is initially (at t = 0) at a given position, namely, x(0) = x 0 , and thus C = x 0. The unique solution, considering the initial position, is:

x(t) =

∫ (^) t

0

v(t)dt + x 0

and x(t) =

v(t)dt + C is the general solution. When both the ODE and the initial condition are given, we say that we have a Cauchy problem or an initial value problem. The physical model of the movement is more meaningful in terms of the l’acceleration a(t): the (second order) ODE representing the movement is:

d^2 x dt^2

= a(t)

Second order equations are usually solved through two steps, each first order. Namely, if we call dx dt(t )= v(t) then the second order equation can be represented as a first order equation in v:

dv(t) dt

= a(t)

Solving this through separation of variables:

v(t) =

∫ (^) t

0

a(t)dt + v 0

and reintroducing dx dt(t )= v(t) one has:

dx(t) dt

= v(t) =

∫ (^) t

0

a(t)dt + v 0

which results through separation of variables in:

x(t) =

∫ (^) t

0

∫ (^) t

0

a(t)dt + v 0 t + x 0

If a(t) were constant, the known expression for the uniformemly accelerated movement results:

x(t) =

at^2 + v 0 t + x 0

The same equations and solutions hold for the more realistic situation of three- dimensional movement, considering that x(t), v(t) and a(t) are 3D points or 3D vectors.

1.5.2 Newton’s Cooling Law

Newton formulated the following Law of Cooling: The rate of loss of heat by a body is directly proportional to the temperature difference between of the body and its surroundings, provided the difference is small.

CHAPTER 1. INTRODUCTION AND BASIC CONCEPTS 7

This law is represented through a DE: t is the time; T (t) the body temper- ature at t; the cooling, represented by C(t), is the rate of change of the body temperature; the initial body temperature is T 0 ; the surroundings temperature is TS. The ODE describing this model is:

C(t) =

dT dt

= c(TS − T (t))

where c is the cooling constant. Through separating the variables T and t one obtains:

T (t) = TS − (TS − T 0 )e−ct

Remark: What happens if TS > T 0?

Stationary solutions

Let us impose dTdt = 0. Then c(TS − T ) = 0 and as c 6 = 0, TS − T = 0 and this equation has the solution T = TS. If a DE has a solution which does not depent on time it is called a stationary solution of the differential equation. These type of solutions make the velocity of change 0, as seen. Thus, if the initial temperature is TS , it will not change, will remain stationary. Any initial temperatute tends to reach the stationary one.

Exercise 3 Using Newton’s cooling law, T ′(t) = k(TS − T (t)):

• The temperature of some food before introducing it in a freezer at − 10 degrees is 25 degrees. After 8 minutes the food is at 0 , when will it reach − 5?
• Which is the stationary solution of the equation?

1.6 Euler’s Method for First Order ODEs

1.6.1 Introduction

The numerical solution of ODEs is also known as numerical integration of ODEs. A basic method of numerical solution we introduce is the finite diferences method. Let us start with a simple case, where the ODE is:

dx dt

= x′(t)

and the initial value is x(0) = x 0. By the definition of derivative at a point:

x′(t) = lim ∆t→ 0

x(t + ∆t) − x(t) ∆t

If we take a small ∆t, one could write, approximately:

x′(t) '

x(t + ∆t) − x(t) ∆t

CHAPTER 1. INTRODUCTION AND BASIC CONCEPTS 9

• that the error (the difference between the exact and the approximate so- lution) is increasing with every step.

Exercise 4 Consider the ODE x′(t) = x(t) with initial value x(0) = 1.

• Check that et^ is the solution of this IVP.
• Compare the values obtained through Euler’s method at times t = 0. 1 , 0. 2 , ..., 1 (i.e., using ∆t = 0. 1 ) to the exact values (et) at those times.
• Compute the absolute and relative errors in the interval t ∈ [0, 1].

1.6.2 Formalisation of Euler’s Method

Consider a slightly more general first order ODE, where the derivative can be isolated, and initial condition:

dx dt = x′(t) = f (t, x) x 0 = x(0)

To apply the numerical method it is convenient to improve the notation. Let us denote ∆t = h, kh = tk, and xk the approximate value of x(tk), the exact solution. Then, the approximation that we wrote above can be expressed as:

xk+1 = xk + hf (tk, xk) x 0 = x(0)

This is an iterative formulation, where the computation of the approximation in the next time requires the approximation in the previous one. This is the (forward) Euler’s method.

Exercise 5 Consider the ODE x′(t) = x^2 (t) + 2t − t^4 with initial value x(0) = 0.

• Check that t^2 is the solution of this IVP.
• Compare the values obtained through Euler’s method at times t = 0. 1 , 0. 2 , ..., 1 (i.e., using ∆t = h = 0. 1 ) to the exact values at those times.
• Determine the errors in the interval t ∈ [0, 1].

1.6.3 Basic Error Concepts

The source of error of Euler’s method is that the exact derivative, f (t, x(t)),

is replaced by an approximation, the fraction x(t+h h)− x(t). Then the local dis- cretization error (of Euler’s method) at each t - thus, the local - as:

L(t, h) = x(t + h) − x(t) h

− f (t, x(t))

Indeed, if we consider only one iteration, the error made, as difference between the exact and the approximate solutions, would be:

x(tk+1) − xk+1 = x(tk+1) − x(tk) − hf (tk, x(tk)) = hL(tk, h)

CHAPTER 1. INTRODUCTION AND BASIC CONCEPTS 10

Thus, the error at a single iteration depends on L(t, h). In fact, x(tk)) is not known, what has been previously computed is its ap- proximation xk; the global discretization error is the error accumulated in n iterations. In a later chapters, errors will be studied in more detail, and more precisely: for instance, showing that under some conditions, the global dis- cretization error is 0(h), i.e., tends to 0 when h tends to zero.

1.7 Supplementary Material

1.7.1 Population Dynamics

Malthus’ model

The Malthus’ model postulates that the population grows proportionally to its current size. If u(t) is the population Malthus’ model can be represented through the ODE:

du dt

= ku

Separating variables one gets:

du u

= kdt

and thus

log u = kt + C u = ekt+C

If the initial population is u(0) = u 0 we obtain

u = u 0 ekt

This model results into an exponential growth of the population.

Logistic growth

The logistic growth model (or Verhulst’s model) is represented by the ODE:

du dt

= ku(1 − u)

and initial condition u(0) = u 0. Separating variables:

du u(1 − u)

= kdt

which can be integrated descomposing the fraction as:

1 u(1 − u)

A

u

B

1 − u

Identifying coefficients, one obtains A = B = 1, which leads to:

log

u 1 − u

= kt + C

u 1 − u

= ekt+C

CHAPTER 1. INTRODUCTION AND BASIC CONCEPTS 12

¨y =

m

(T −

cρsv^2 )sinθ − m˙ m

y˙ − g

where: x˙ = dxdt , ¨x = d

(^2) x dt^2 , ˙y^ =^

dy dt , ¨y^ =^

d^2 y dt^2 ,^ m˙^ =^

dm dt , ¨m^ =^

d^2 m dt^2. Indeed it is a system of coupled second order DEs.

Exercise 6 Find the general solution of the ODE x′−t = sin(t) using separation of variables.

Exercise 7 The dynamics of an infectious bacterial population P (t) within a body evolves according to the ODE:

P ′(t) = (kcos(t))P (t), k > 0

• Obtain a general solution of the population evolution.
• If k = 1 and the initial population is 100, provide the dynamics of the population.
• Sketch the graphical representation of this solution.

Exercise 8 Attempt to find the general solutions of some of the ODEs in exercise 1 through separation of variables.

Chapter 2

State, Simulation,

Visualization, Processing

We shall use the programming environment Processing in this course, because it makes the visualization especially simple. A physical model is represented through DEs involving several variables: those representative of the model are called the state representation. Through numerical methods, the DEs can be approximately solved, i.e., the system can be simulated. The simulation can be visualised: in this chapter we introduce the state visualisation. Simulation is discussed in depth in another chapter.

2.1 State

2.1.1 Representing the State

The first step towards the visualization of a (physical) system is to decide the state description, i.e., the variables representing the system at each time, in a way both complete - without missing anything needed - and economical - withouth unnecessary data, which would need extra computation at each time -.

2.1.2 A First Example: the Pendulum

Assume we wish to simulate the pendulum movement restricted to two dimen- sions. Two variables are sufficient, namely, its angle with respect to the equi- librium at each time, θ; and the angular speed θ′. An extra fixed magnitude is needed: the pendulum arm length. The following image, taken from Wikipedia, shows different variables, some of them are enough to represent the state..

CHAPTER 2. STATE, SIMULATION, VISUALIZATION, PROCESSING 15

float DX = 0.1;

// Size of the array of heights representing the surface wave in this // 10 m. world each time. int ArraySize = ( int )ceil( WorldSize / DX );

// Definition of the heights array. float[] StateHeightArray = new float[ArraySize];

This is not the complete state necessary for the water surface simulation. It will be completed later when discussing the wave equation (a PDE). Remark: the array dimension is obtained through a formula, which is more flexible: the deposit size and the horizontal sampling (resolution) - in this case 0.1 m. An alternative formulation, where the resolution is obtained from world size and array dimension would be:

// World size, in metres. float WorldSize = 10.0;

// Array size int ArraySize = 100;

// Separation between two consecutive height measures, in metres. float DX = WorldSize / ( float )ArraySize;

// Definition of the heights array. float[] StateHeightArray = new float[ArraySize];

Before discussing visualization, let us introduce basic programming with Pro- cessing.

2.2 Basic Processing

Processing was created to enable quick simulation and visualization, que esta pensat per simular i visualitzar molt rapidament. What follows is based on: Reas, C, and Fry, B.: Getting Started with Processing. You should download and install Processing first.

2.2.1 Basic Visualization

Example 1

Open Processing and try:

// Draw a circle centre (50,50) and radius 80. ellipse(50, 50, 80, 80);

You should have copied the code and pasted it in the editor. Press Run to visualize the results. Try next:

ellipse(30, 40, 20, 30);

CHAPTER 2. STATE, SIMULATION, VISUALIZATION, PROCESSING 16

and

ellipse(30, 40, 30, 20);

Remark:

• A circle is drawn as an ellipse with two equal axes.
• The origin (of the pixel coordinate sysem used) is at the top left of the window. Positive values of the first coordinate go to the right, those of the second, downwards.
• A semicolon marks the end of a line; while comments (beginning with //) are ignored, and are useful to document the code.

Example 2

The next example is interactive:

void setup() { // Instead of the default window size, use a custom window size. // Our (rectangular) window is 480 pixels width and 120 height size(480, 120); smooth(); }

void draw() { if (mousePressed) { // If the mouse button is pressed black (0) is painted fill(0); } else { // If the mouse button is not pressed white (255) is painted fill(255); } // What is being painted (either black or white) is a circle centered // where the mouse lies and 80 pixel radius ellipse(mouseX, mouseY, 80, 80); }

The mini-programs are called sketches). Open a saved sketch, create a new one, export a sketch.

Drawing Elementary Geometric Shapes

A line segment is defined by its two endpoints; a triangle by its 3 vertices, a rectangle by 4 (thus, 8 valus). Note that a square is a specific case of a rectangle. Essay, for instance:

line(20, 50, 420, 110); triangle(347, 54, 392, 9, 392, 66); triangle(158, 55, 290, 91, 290, 112); quad(158, 55, 199, 14, 392, 66, 351, 107);

CHAPTER 2. STATE, SIMULATION, VISUALIZATION, PROCESSING 18

• Using colour:

size(480, 120); noStroke(); smooth(); background(0, 26, 51); fill(255, 0, 0); // Red circle ellipse(132, 82, 200, 200); fill(0, 255, 0); // Green circle ellipse(228, -16, 200, 200); fill(0, 0, 255); // Blue circle ellipse(268, 118, 200, 200);

By default the most recently painted circle covers the previous one.

2.2.2 Variables (and their Visualization)

Basic Aspects

To declare and assign values to variables

// x is declared as an integer variable int x; // Assign a value to x x = 12;

More concisely

// x is declared as an integer variable and a value is assigned int x = 12;

There are default variables, which can be reused. The following code paints circles and lines depending on the window size:

size(480, 120); smooth(); // segment from (0,0) to (480, 120) line(0, 0, width, height); // segment from (480, 0) to (0, 120) line(width, 0, 0, height); // circle centred at (240, 60) of radius 30 ellipse(width/2, height/2, 30, 30);

We can perform operations (+,-,*,/) with the variables: see the example

size(480, 120); int x = 25; int h = 20; int y = 25; // upper

CHAPTER 2. STATE, SIMULATION, VISUALIZATION, PROCESSING 19

rect(x, y, 300, h); x = x + 100; // middle rect(x, y + h, 300, h); x = x - 250; // lower rect(x, y + h*2, 300, h);

Declaration and assignment can use operations:

// Assigning 24 to x int x = 4 + 4 * 5;

Remark that in the state definition some variables were declared as integer int, or as real float, and operations were used too.

Loops (for)

One can use loops for repetitive tasks, for instance the following code, which paints lines:

size(480, 120); smooth(); strokeWeight(8); line(20, 40, 80, 80); line(80, 40, 140, 80); line(140, 40, 200, 80); line(200, 40, 260, 80); line(260, 40, 320, 80); line(320, 40, 380, 80); line(380, 40, 440, 80);

can be replaced by

size(480, 120); smooth(); strokeWeight(8); for (int i = 20; i < 400; i += 60) { line(i, 40, i + 60, 80); }

Note assignment i = i +60, shortened to i += 60. Within the signs { } we have a block; after for we initialize, then test the condition and update; while the test gives true, a statement is executed, and if false, we exit the loop.

2.2.3 Two Basic Methods: Setup and Draw

It is easy to see how both methods work through a simple example, using a function that writes on the console, println():