## Busca en el resumen del documento

Numerical Methods

Josep Blat

Universitat Pompeu Fabra

Contents

1 Introduction and Basic Concepts 3 1.1 Basic Concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Orientation of the Course . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Basic ODE Concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Analytical Solutions through Separation of Variables . . . . . . . . . . . . . . . 5 1.5 Models and Corresponding ODEs . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.5.1 Trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.5.2 Newton’s Cooling Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.6 Euler’s Method for First Order ODEs . . . . . . . . . . . . . . . . . . . . . . . 7 1.6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.6.2 Formalisation of Euler’s Method . . . . . . . . . . . . . . . . . . . . . . 9 1.6.3 Basic Error Concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

1.7 Supplementary Material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.7.1 Population Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.7.2 The Trajectory of a Rocket . . . . . . . . . . . . . . . . . . . . . . . . . 11

2 State, Simulation, Visualization, Processing 13 2.1 State . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.1.1 Representing the State . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.1.2 A First Example: the Pendulum . . . . . . . . . . . . . . . . . . . . . . 13 2.1.3 Second Example: Water Surface Simulation . . . . . . . . . . . . . . . . 14

2.2 Basic Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2.1 Basic Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2.2 Variables (and their Visualization) . . . . . . . . . . . . . . . . . . . . . 18 2.2.3 Two Basic Methods: Setup and Draw . . . . . . . . . . . . . . . . . . . 19

2.3 Visualizing Pendulum and Waves . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.1 Pendulum Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.2 Drawing Surface Waves . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.4 Other Examples, Coding Style . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.4.1 Coding Conventions Used . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.2 Abstracting the State . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3 The Elastic Spring: One-step Methods 25 3.1 A Physical Model of the Spring . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Numerical Solution of the Spring ODE through Euler . . . . . . . . . . . . . . 26

3.2.1 From a Second Order ODE to a First Order System . . . . . . . . . . . 26 3.2.2 Euler’s Method Applied to the First Order System . . . . . . . . . . . . 26

3.3 Coding Euler’s Method Applied to the 1D Spring . . . . . . . . . . . . . . . . . 26 3.3.1 State and Visualization of the 1D Spring . . . . . . . . . . . . . . . . . . 26 3.3.2 The Time Step Function . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3.3 Comparing to the Exact Solution . . . . . . . . . . . . . . . . . . . . . . 29

3.4 Higher Order, Implicit, Multistep Methods . . . . . . . . . . . . . . . . . . . . 31 3.5 A Semi-implicit Method on the Spring: Euler-Cromer . . . . . . . . . . . . . . 31 3.6 Higher Order Methods: Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . 32

3.6.1 Runge-Kutta 2 (Heun) Method . . . . . . . . . . . . . . . . . . . . . . . 32 3.6.2 Runge-Kutta 4 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.7 Simulation of the Elastic Spring with RK2 . . . . . . . . . . . . . . . . . . . . . 34 3.8 Simulació de la molla utilitzant el mètode Runge-Kutta d’ordre 4 (RK4) . . . . 35

1

CONTENTS 2

3.9 Supplementary Material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.9.1 Analytical Solution of the Spring ODE . . . . . . . . . . . . . . . . . . . 36 3.9.2 Variants of the Spring Model . . . . . . . . . . . . . . . . . . . . . . . . 38

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.

3

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.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 x0 = x(0) exists and is unique, for t ∈ (0, a). A successive

CHAPTER 1. INTRODUCTION AND BASIC CONCEPTS 5

approximation technique is used to prove this theorem. When the initial value is not fixed, the so called general solutions contain

one or more arbitrary constants.

Exercise 1 Couple each ODE with its general solution in the following table, using the definition of solution. Remark the relationship between order of the DE and number of arbitrary constants.

1 x′ = 0 A x = t2 + t+ C 2 x′ = 2t+ 1 B x = 2t2 + C1t+ C2 3 x′ = x C x = −5 + Ce−2t 4 x′ + 2tx = t D x = C1sin(2t) + C2cos(2t) 5 x′ + 2x+ 10 = 0 E x = C

6 x′′ + 4x = 0 F x = 12 + Ce −t2

7 x′′ = 4 G x = Cet

1.4 Analytical Solutions through Separation of Variables

If the DE can be transformed to the form f(x)dx = g(t)dt, i.e., the variables are separated, the analytical general solution of the equation is

∫ f(x)dx =∫

g(t)dt+ C.

Exercise 2 Find the general solution of the ODE x′ = x using separation of variables.

1.5 Models and Corresponding ODEs

1.5.1 Trajectories

The movement of a particle along a line, the X axis for instance, can be rep- resented through a first order ODE, where v(t) is the velocity or speed at each instant t:

dx(t)

dt = v(t)

The analytical general solution of this equation can be obtained through vari- ables separation. It can be rewritten as:

dx = v(t)dt

and integrating both sides of the equal sign:

x(t) =

∫ v(t)dt+ C

If v(t) were constant, x(t) = vt+C; depending on the complexity of v(t) a closed solution might not exist. It seems that there are infinitely many solutions, each for a value taken by C. But considering the physical model, the solution is in

CHAPTER 1. INTRODUCTION AND BASIC CONCEPTS 6

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

x(t) =

∫ t 0

v(t)dt+ x0

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:

d2x

dt2 = a(t)

Second order equations are usually solved through two steps, each first order.

Namely, if we call dx(t)dt = 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+ v0

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

dx(t)

dt = v(t) =

∫ t 0

a(t)dt+ v0

which results through separation of variables in:

x(t) =

∫ t 0

∫ t 0

a(t)dt+ v0t+ x0

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

x(t) = 1

2 at2 + v0t+ x0

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 T0; 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 − T0)e−ct

Remark: What happens if TS > T0?

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) = x0. 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 8

If we replace ' by = (an approximation is made, and an error is introduced):

x′(t)∆t = x(t+ ∆t)− x(t)

x(t+ ∆t) = x(t) + x′(t)∆t

The previous formula, with a fixed small ∆t can be used as a method to compute (approximate) values of the solution at instants n∆t:

• start with the initial value x(0) = x0

• compute the next value

x(∆t) = x(0) + ∆tx′(0)

• and iteratively:

x(n∆t) = x((n− 1)∆t+ ∆t) = x((n− 1)∆t) + x′((n− 1)∆t

The following image (taken from http://tutorial.math.lamar.edu/) visualises the computation of the first step, showing the approximate solution, and the exact solution, where the unknown is called y instead of x.

The following image (taken from the Wikipedia) sketches the computation of four steps, showing the approximate solution in red - interpolating linearly between the four computed values -, and the exact solution in blue.

Remark that

• the exact and approximate solutions are only equal at the initial time,

• the approximation used to compute the next value is based on the tangent

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) x0 = 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) x0 = 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) = x2(t) + 2t− t4 with initial value x(0) = 0.

• Check that t2 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)−x(t)h . 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) = u0 we obtain

u = u0e kt

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) = u0. 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 11

and thus:

u = ekt+C(1− u) u = e kt+C

1 + ekt+C

Visualising u versus t a characteristic sigmoidal curve appears (the logistic curve), instead of an exponential growth.

If we impose dudt = 0, then u(1− u) = 0 and this equation has two solutions, u = 0 and u = 1, which are the stationary solutions of the differential equation. These solutions make the velocity of change 0, as seen. Thus, if the initial population is 0 or 1, the population will not change, will remain stationary.

Any nonzero initial population tends to reach the population limit u = 1. Indeed, if the initial population is larger than 1, it decreases; if it is smaller, it grows.

Predator-prey models

The Lotka-Volterra model considers two populations: the prey u(t) and the predator v(t); the prey grows proportionally to its size, but decreases propor- tionally to the size of the predator and its own size; the predator decreases proportionally to its size, and grows proportionally to the size of the prey and its own size. This can be represented through the system of DEs:

du

dt = αu+ βuv

dv

dt = γv + δuv

where α > 0, β < 0, γ < 0 and δ > 0 (positive sign means growth, negative sign decay).

1.7.2 The Trajectory of a Rocket

A rocket can fly long distances due to the consumption of its own fuel. This means that its mass varies during the flight. A relatively simple model of the trajectory of a rocket takes into account only three forces: a vertical one due to gravity ; the thrust of the rocket T assumed to be constant while the fuel lasts; and friction of the atmosphere, proportional to the rocket speed v and its section s but negative. It does not consider the thrift due to the wind, for instance. Other simplifications are to assume that the trajectory is not too long nor too high, so that the Earth can be considered as plane, the air density and the gravity as constant.

Then, the trajectory is planar, so that the rocket position can be represented by means of just two coordinates, x horizontal and y vertical; similarly velocity through vx and vy, acceleration ax and ay acceleració vertical; θ is the angle between ~v(t) and a horizontal line; and ρ is the air density. Using the law of mechanics

d(m~v)

dt = ~F

the following ODEs are obtained:

ẍ = 1

m (T − 1

2 cρsv2)cosθ − ṁ

m ẋ

CHAPTER 1. INTRODUCTION AND BASIC CONCEPTS 12

ÿ = 1

m (T − 1

2 cρsv2)sinθ − ṁ

m ẏ − g

where: ẋ = dxdt , ẍ = d2x dt2 , ẏ =

dy dt , ÿ =

d2y dt2 , ṁ =

dm dt , m̈ =

d2m dt2 . 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..

13

CHAPTER 2. STATE, SIMULATION, VISUALIZATION, PROCESSING 14

The following Processing code can represent the pendulum state:

// Length of the pendulum arm, in metres. In fact, it is a constant

// which does not belong really to the de simulation state

float PendulumLength;

// Angle of the pendulum at each instant, in radians

float StatePendulumTheta;

// Variation of the angle in time, angular velocity

// at a specific time, in radians per second

float StatePendulumDthetaDt;

Remark: the magnitudes should be represented through appropriate physical measure units, as in the examples.

2.1.3 Second Example: Water Surface Simulation

Another more complex example, which we shall study in detail this course, is related to the simulation of waves on the water surface. To simplify, consider the 2D section of a water deposit (see the next image). The height is the main simulation parameter to simulate and visualize.

The vertical height corresponding to each horizontal position at each time is what should be simulated/visualized. In practice one cannot represent infinite horizontal positions; if, for example, the horizontal dimension measures 10 m we could represent the heightes every 10 (horizontal) cm. The state in Processing:

// World size, in metres.

float WorldSize = 10.0;

// Space between two consecutive heights, in metres (.1 m = 10 cm).

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 està pensat per simular i visualitzar molt ràpidament. 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 17

Parts of an ellipse can be drawn through giving the start and end angles (in radians!). Trying the following examples illustrates the positive sign as well:

arc(90, 60, 80, 80, 0, HALF_PI);

arc(190, 60, 80, 80, 0, PI+HALF_PI);

arc(290, 60, 80, 80, PI, TWO_PI+HALF_PI);

arc(390, 60, 80, 80, QUARTER_PI, PI+QUARTER_PI);

The function radians(-) transforms degrees into radians (example: radians(90) = π 2 ).

Less Basic Drawing

Find a tutorial with examples to:

• control the drawing order (paint a rectangle on an ellipse and viceversa)

• smooth lines (function smooth( ), already used; noSmooth( ) will disable smooth)

• define the thickness of the stroke (strokeWeight( );) and more advanced functions

size(480, 120);

smooth();

ellipse(75, 60, 90, 90);

// stroke thickness 8 pixels

strokeWeight(8);

ellipse(175, 60, 90, 90);

ellipse(279, 60, 90, 90);

// stroke thickness 20 pixels

strokeWeight(20);

ellipse(389, 60, 90, 90);

• define colour (grey level or RGB)

size(480, 120);

smooth();

// Black

background(0);

// Light grey

fill(204);

// Circle in light grey

ellipse(132, 82, 200, 200);

// Middle grey

fill(153);

// Circle in middle grey

ellipse(228, -16, 200, 200);

// Dark grey

fill(102);

// Circle in dark grey

ellipse(268, 118, 200, 200);

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():

CHAPTER 2. STATE, SIMULATION, VISUALIZATION, PROCESSING 20

void setup()

{

println("I’m starting");

}

void draw()

{

println("I’m running");

}

Alternating quickly Run and Stop, one sees that I’m starting is written once, while I’m running is written until we stop. Within setup one usually puts things that are defined just once such as the window size definition; within draw things such as the user’s input, as the program should normally watch continuously for it.

The second example given can be better understood now: at set up the window was defined and the fixed smooth function; at draw circles are painted, centered where the user had the mouse and the colour determined by the mouse button being pressed or not. The example with comments follows:

void setup()

{

// Definition of the (rectangular) window size. Width 480 pixels, Height 120.

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);

}

2.3 Visualizing Pendulum and Waves

2.3.1 Pendulum Visualization

First, we transform the real coordinate systema to the screen system, which uses pixels. If a 2 m world is visualized on a 300x300 screen, then the heading of the program could be:

int WindowWidthHeight = 300;

float WorldSize = 2.0;

float PixelsPerMeter;

CHAPTER 2. STATE, SIMULATION, VISUALIZATION, PROCESSING 21

float OriginPixelsX;

float OriginPixelsY;

And within setup the transformation values are computed. Within the X origin is set at the window middle, while the Y origin relatively high, in the upper quadrant:

PixelsPerMeter = (( float )WindowWidthHeight ) / WorldSize;

OriginPixelsX = 0.5 * ( float )WindowWidthHeight;

OriginPixelsY = 0.25 * ( float )WindowWidthHeight;

Within draw we call a function DrawState to draw some of the state ele- ments: the arm, the pendulum pivot and bob; at each time we need the angle as well. The code follows:

// The DrawState function assumes that the simulation coordinate system

// is measured and the pendulum pivot is at the origin. Draw arm, pivot, bob

void DrawState()

{

// Compute the arm end.

float ArmEndX = PixelsPerMeter * PendulumLength * sin( StatePendulumTheta );

float ArmEndY = PixelsPerMeter * PendulumLength * cos( StatePendulumTheta );

// Draw the arm.

strokeWeight( 1.0 );

line( 0.0, 0.0, ArmEndX, ArmEndY );

// Draw the pivot.

fill( 0.0 );

ellipse( 0.0, 0.0,

PixelsPerMeter * 0.03,

PixelsPerMeter * 0.03 );

// Draw the bob.

fill( 1.0, 0.0, 0.0 );

ellipse( ArmEndX, ArmEndY,

PixelsPerMeter * 0.1,

PixelsPerMeter * 0.1 );

}

The following image shows the visualization assuming the angle is π/6, if the program were correct. This way one can test whether the visualization code is working. In order to simulate, we need to introduce the solution of the corresponding DE.

CHAPTER 2. STATE, SIMULATION, VISUALIZATION, PROCESSING 22

2.3.2 Drawing Surface Waves

The next example, simulated later, is the visualization of the heights 1D field for the surface waves. Now we need a more complex DrawState function which paints blue rectangles whose width is the resolution, and whose height is that of the current state; they should be painted one next to each other filling the whole deposit. At the beginning the suitable world to screen coordinate system transformations are defined. The code can be:

void DrawState()

{

float OffsetY = 0.5 * ( float )WindowHeight;

for ( int i = 0; i < ArraySize; ++i )

{

float SimX = DX * ( float )i;

float PixelsX = ( float )( i * PixelsPerCell );

float SimY = StateHeight[i];

float PixelsY = SimY * (( float )PixelsPerCell ) / DX;

float PixelsMinY = OffsetY - PixelsY;

float PixelsHeight = (( float )WindowHeight) - PixelsMinY;

fill( 0.0, 0.0, 1.0 );

rect( PixelsX, OffsetY - PixelsY, PixelsPerCell, PixelsHeight );

}

}

The visualization, if correctly programmed would be as in the following image. To make it move, we need the simulation, the solution of the DE.

CHAPTER 2. STATE, SIMULATION, VISUALIZATION, PROCESSING 23

2.4 Other Examples, Coding Style

More complex states can be illustrated through the following example. There are a number of particles, with several attributes each, such as position, mass, speed, etc. The state at a time depends usually on the values at previous times, and thus the latter should be part of the state, to facilitate the simulation. The following code exemplifies this:

// Max number of particles.

int MAX_NUM_PARTICLES = 1024; int ArraySize = MAX_NUM_PARTICLES;

// Number of active particles.

int N;

// Arrays of data of each particle. There is an array for each attribute,

// separating eventually the vector components.

// This illustrates a simulation strategy:

// structs of arrays, rather than arrays of structs.

float[] StateMass = new float[ArraySize];

float[] StateRadius = new float[ArraySize];

float[] StatePosX = new float[ArraySize];

float[] StatePosY = new float[ArraySize];

float[] StateVelX = new float[ArraySize];

float[] StateVelY = new float[ArraySize];

float[] StatePrevPosX = new float[ArraySize];

float[] StatePrevPosY = new float[ArraySize];

The following could be the full state of the surface waves example:

// World size, in metres.

float WorldSize = 10.0;

// Distància entre dues mesures consecutives d’altura en metres.

// (.1 metres equival a 10 cm)

float DX = 0.1;

CHAPTER 2. STATE, SIMULATION, VISUALIZATION, PROCESSING 24

// Array size of the heights representing the surface wave

// of the 10 m. world at some time.

int ArraySize = ( int )ceil( WorldSize / DX );

// Defining the heights array.

float[] StateHeight = new float[ArraySize];

float[] StateDheightDt = new float[ArraySize];

float[] StatePressure = new float[ArraySize];

float[] StateFloorHeight = new float[ArraySize];

float[] StatePrevHeight = new float[ArraySize];

2.4.1 Coding Conventions Used

Coding conventions are important to facilitate use and re-use of code. Remark that in each example so far always:

• the size of arrays is ArraySize

• the names of the state variables have the prefix State,

• the size of the world is WorldSize,

• the horizontal resolution is DX,

and so on.

2.4.2 Abstracting the State

Another abstraction level can come from putting all the magnitudes within a state at each time; an index will denote each magnitude used. Thus each magnitude could be directly accessed if required. This macro-vector will be always called State Vector. See the particles example:

int ArraySize = 1024;

int NumStateArrays = 8;

float[][] State = new float[NumStateArrays][ArraySize];

int StateIndexMass = 0;

int StateIndexRadius = 1;

int StateIndexPosX = 2;

int StateIndexPosY = 3;

int StateIndexVelX = 4;

int StateIndexVelY = 5;

int StateIndexPrevPosX = 6;

int StateIndexPrevPosY = 7;

// To access the position x of particle i:

float xPos = State[StateIndexPosX][i];