In 1925 Umberto d’Ancona determines the percentage of sharks of total amount of fish caught in the harbor of Trieste:

(1)

1 Ordinary Differential Equations: Examples, theoretical foundations and tools

1.1 Volterra’s principle - a motivating example

In 1925 Umberto d’Ancona determines the percentage of sharks of total amount of fish caught in the harbor of Trieste:

1914 1915 1916 1917 1918 1919 1920 1921 1922 1923 0

5 10 15 20 25 30 35 40

Remarkable: reduced fishing (World War I) a disadvantage for food fish?

ODE basics TU Bergakademie Freiberg

6 The situation can be described by a Predator-Prey model (Vito Volterra, 1860-1940):

x(t) . . . prey population at time t (food fish), y(t) . . . predator population at time t (sharks).

In the absence of predators the prey population would increase according to the Malthusian

^∗

law of growth

x

⁰

(t) = a x(t) (with a constant a > 0), i.e., increase proportional to current number (exponential growth)

x(t) = x(0) e ^at for t ≥ 0

(applies if population not too dense and sufficient food available).

Thomas Robert Malthus (1766-1834)

7 Decrease of prey and growth of predators are proportional to the number of predator-prey contacts (per unit time), i. e. to

x(t)y(t).

Hence, a mathematical model (without fishing) is given by the following system of (ordinary) differential equations.

x

⁰

(t) = a x(t) − b x(t)y(t) (a, b > 0),

y

⁰

(t) = −c y(t) + d x(t)y(t) (c, d > 0). (1) For uniqueness of the solution we have to set an initial condition, e. g.

measured values for x and y at time t = 0.

(2)

It can be shown that the solutions are periodic, i. e.

x(t) = x(t + T ), y(t) = y(t + T ),

for some T > 0, and that the means of population numbers fulfill

¯ x = 1

T Z T

0 x(t) dt = c

d , und y ¯ = 1 T

Z T 0

y(t) dt = a b . For the solution of the ODE one usually applies numerical algorithms.

Typical solutions look as follows:

a/b c/d

x(t)

y(t)

9 To take account of moderate fishing (0 < e < a), one modifies the system (1) to

x

⁰

(t) = a x(t) − b x(t)y(t) − e x(t) = (a − e) x(t) − b x(t)y(t), y

⁰

(t) = − c y(t) + d x(t)y(t) − e y(t) = − (c + e) y(t) + d x(t)y(t). (2) Structurally it’s the same system as (1), but with coefficients a − e instead of a and c + e instead of c.

For the new means we get

¯

x F = c + e d > c

d = ¯ x und y ¯ F = a − e b < a

b = ¯ y.

Volterra’s principle

Moderate fishing (e < a) increases the average number of food-fish and reduces the average number of sharks.

10 A phase plot of the solution, i. e. the plot of the curve t 7→

x(t) y(t)

, 0 ≤ t ≤ T simplifies the comparison of both solutions:

c/d (c+e)/d (aïe)/b

a/b

Prey

Predator

without fishing with fishing

(3)

General context

The equations (1) form a system of ordinary differential equations, i. e.

variation in time, not in space;

at each point in time system state characterized by finitely many numbers.

Of course we need precise definitions!

Problems leading to ODEs are at least as old as Newton (celestial mechanics) and Euler (brachistochrone).

But solutions for ODEs from real problems typically are only possible with numerical methods.

12 1.2 Theoretical Background An expression of the form

F

t, y, y

⁰

, y

⁰⁰

, . . . , y ⁽ⁿ⁾

= 0 (3)

with a function F : R ⁿ⁺² ⊃ M → R is called an ordinary differential equation (ODE) of order n.

A function y : R ⊃ I → R is called solution von (3) on the interval I if it is n times continuously differentiable on I (y ∈ C ⁿ (I)) and

F

t, y(t), y

⁰

(t), y

⁰⁰

(t), . . . , y ⁽ⁿ⁾ (t)

= 0 for all t ∈ I.

13 Comments on terminology

The ODE (3) has order n, because n is the order of the highest occurring derivative.

The term ordinary refers to the fact that all derivatives of y are taken with respect to the same variable t.

The ODE (3) is called implicit, because of the zero on the right hand side.

Opposed to this we can consider an explicit ODE of order n, where the derivative of highest order is isolated on the left hand side:

y ⁽ⁿ⁾ = f

t, y, y

⁰

, . . . , y ⁽ⁿ⁻¹⁾

.

(4)

We shall consider exclusively systems of explicit first-order ODEs.

(We shall see shortly why first-order systems are general enough) y ₁

⁰

= f 1 (t, y 1 , y 2 , . . . , y n )

y ₂

⁰

= f 2 (t, y 1 , y 2 , . . . , y n ) ... = ...

y

⁰

_n = f n (t, y 1 , y 2 , . . . , y n )

(DE)

with the n unknown functions y 1 , y 2 , . . . , y n .

Any system of n functions y 1 = y 1 (t), . . . , y n = y n (t) ∈ C ¹ (I) which satisfies (DE) for all t ∈ I is called a solution of (DE) on I.

15 Example. The system y ₁

⁰

= 1 y ₂

⁰

= 2y 1

possesses the solutions

y 1 (t) = t + α,

y 2 (t) = t ² + 2αt + β (α, β ∈ R ) on (−∞, ∞).

Uniqueness of solution requires specifying initial conditions, e. g.

y 1 (0) = 1, y 2 (0) = 2.

With these additional conditions y 1 (t) = t + 1, y 2 (t) = t ² + 2t + 2, is the unique solution.

16 In general: The problem of finding a solution of (DE) which also satisfies the initial condition(s)

y 1 (t 0 ) = y 0,1 , . . . , y n (t 0 ) = y 0,n (AB) is known as an initial value problem (IVP) for the ODE (DE).

Setting

y :=



  y 1

...

y n



  , f :=



  f 1

...

f n



  , y 0 :=



  y 0,1

...

y 0,n



  ,

we obtain the short notation

y

⁰

= f (t, y), (4a)

y(t 0 ) = y 0 . (4b)

(5)

Transformation of higher order ODEs

ODEs of higher order may be transformed into an equivalent system of first-order ODEs: as an example, y

⁰⁰⁰

+ 3y

⁰⁰

+ y

⁰

= sin(t) becomes



 y

⁰

₁ y

⁰

₂ y

⁰

₃



 =



 y 2

y 3

−3y 3 − y 2 + sin(t)



 (y 1 = y), (y 2 = y

⁰

₁ = y

⁰

), (y 3 = y

⁰

₂ = y

⁰⁰

).

The following IVP is a model for a one dimensional damped oscillator (Stokes type frictional force):

m¨ x = − kx − δ x, ˙ x(0) = x 0 , x(0) = 0 ˙ (5) (m, k, δ are positive constants). Numerical methods for ODEs are usually formulated for first-order ODE systems. Make (5) “ready to use“ by performing a transformation as above.

18 The explicit dependence of the right-hand side on the independent variable (here t) may be removed by adding one further component

y 4 (t) = t (d.h. y

⁰

₄ (t) = 1), y 4 (t 0 ) = t 0

to the solution vector y. This results in an autonomous differential equation y

⁰

= f(y).

In the example above we get

y

⁰

(t) = f (y(t)), f(y) =



 



y 2

y 3

− 3y 3 − y 2 + sin(y 4 ) 1



 

 .

19 Existence and uniqueness of solutions Theorem 1.1 (Picard-Lindelöf

^∗

).

For the IVP

y

⁰

= f (t, y), y(t 0 ) = y 0 , (IVP) assume, that f is continuous in the ‘box’

Q := { (t, y) : | t − t 0 | ≤ a, k y − y 0 k ≤ b } ⊂ R ⁿ⁺¹ . Moreover, let f satisfy the Lipschitz-condition

kf(t, y) − f(t, ˜ y)k ≤ Lky − ˜ yk ∀ (t, y), (t, y) ˜ ∈ Q. (Lip) Then (IVP) has a unique solution on I := [t 0 − α, t 0 + α],

where α := min { a, b/M } and M := max {k f(t, y) k : (t, y) ∈ Q }.

∗

Emile Picard (1856–1941); Ernst Leonard Lindelöf (1870–1946);

Rudolf Otto Sigismund Lipschitz (1832–1903)

(6)

Remarks

(a) For the remainder of the course we assume the fundamental condition (Lip) to hold.

(b) (IVP) possesses a unique solution in [t 0 − a, t 0 + a] if f satisfies the condition (Lip) in

Q ˜ := {(t, y) : |t − t 0 | ≤ a, kyk < ∞}.

(c) If f is continuously differentiable on Q w.r.t. y and f

y

= [∂f i /∂y j ] _1≤,i,j≤n

denotes its Jacobian matrix, then the mean value theorem implies that the requirements of Theorem 1.1 hold with

L := sup

(t,y)

∈

Q k f

y

(t, y) k < ∞ .

21 (d) (IVP) still possesses solutions if one requires only that f be continuous on Q (Existence theorem of Peano

^∗

). These are, in general, no longer unique.

Example: For the IVP y

⁰

= f (t, y) = √

y, y(0) = 0, Q = R × [0, ∞), all functions

y λ (t) =

( 0, 0 ≤ t ≤ λ, (t − λ) ² /4, t ≥ λ, with λ ≥ 0 are solutions.

∗

Giuseppe Peano (1858–1932)

22 Show that the IVP

y

⁰

= λy, y(0) = y 0

has a unique solution on R, no matter how the constants λ, y 0 are chosen. Do you remember this solution from your basic course in higher mathematics?

Of course, the Picard-Lindelöf theorem can not applied to the example

from the last page (otherwise the solution would be unique). What

assumptions are violated?

(7)

Continuous dependence on data Theorem 1.2.

Assume the IVP

y

⁰

= f (t, y), y(t 0 ) = y 0

satisfies the conditions of Theorem 1.1. Assume further that, for a second IVP

z

⁰

= g(t, z), z(t 0 ) = z 0 , its right-hand side g is continuous in Q and there holds

ky 0 − z 0 k ≤ γ and kf (t, y) − g(t, y)k ≤ δ ∀(t, y) ∈ Q.

Then solutions y and z of these IVPs on an interval I satisfy ky(t) − z(t)k ≤ γ e ^L(t−t

⁰

⁾ + δ

L

e ^L(t−t

⁰

⁾ − 1

, ∀t ∈ I.

24 Theorem 1.2 tells you, that small perturbations of the initial condition or of the right hand side will cause only small perturbations in the solution, provided that the elapsed time is small.

This is a very important fact, since in numerical algorithms errors are unavoidable!

Give an upper bound to the relative error in the solution at t = 1, when you consider the IVP

y

⁰

= y + 0.0001 sin t, y(0) = 0.9999 instead of

y

⁰

= y, y(0) = 1.

25 1.3 Matrix Functions

In this section A ∈ C ⁿ

^×

ⁿ shall always denote a square matrix.

For a function f : C ⊇ M → C we explain when and how the matrix f(A) ∈ C ^n×n is defined and review its most important properties.

For ODEs the matrix exponential exp(A), i. e. the exponential function applied to A, is of particular importance.

For some elementary functions f, the definition of the matrix f(A) is obvious. For example, if f ∈ P m is a polynomial of degree m,

f(λ) = α 0 + α 1 λ + α 2 λ ² + · · · + α m λ ^m with coefficients α j ∈ C , j = 0, 1, . . . , m, then

f(A) = α 0 I + α 1 A + α 2 A ² + · · · + α m A ^m ∈ C ⁿ

^×

ⁿ .

(8)

Lemma 1.3 (Properties of p(A) for polynomials p).

Let p be a polynomial of grade m, i. e. p ∈ P m .

(1) If A = diag(A 1,1 , A 2,2 , . . . , A k,k ) is block diagonal with square diagonal blocks A j,j ∈ C ⁿ

^j^×

ⁿ

^j

(j = 1, 2, . . . , k),

n 1 + n 2 + · · · + n k = n, then

p(A) = diag(p(A 1,1 ), p(A 2,2 ), . . . , p(A k,k )).

(2) If T ∈ C ⁿ

^×

ⁿ is invertible and B := T AT

⁻

¹ , then there holds p(B) = T p(A) T

⁻

¹ .

(3) λ is an eigenvalue of A with associated eigenvector v if, and only if, p(λ) is an eigenvalue of p(A) with associated eigenvector v, i.e.

Av = λv ⇐⇒ p(A)v = p(λ)v.

27 Example. We determine m k (J) for the monomial m k (λ) = λ ^k and a Jordan

^∗

block

J = J(λ) =



 

  λ 1

λ 1

... ...

λ 1

λ



 

 

∈ C ^n×n .

Direct calculation shows that m k (J) = J ^k is an upper triangular matrix with Toeplitz structure. The value of the j-th diagonal is

k(k − 1) · · · (k − j + 1)λ ^k

⁻

^j

j! = 1

j!

d ^j

dλ m k (λ) j = 0, 1, . . . , n − 1.

∗

Marie Ennemond Camille Jordan (1838–1922)

28 In other words

m k (J) = J ^k =



 

 

m k (λ) m

⁰

_k (λ) · · · ^m

(n−2)

k

(λ)

(n

−

2)!

m

⁽ⁿ_k⁻¹⁾

(λ) (n

−

1)!

m k (λ) · · · ^m

(n−3)

k

(λ)

(n

−

3)!

m

⁽ⁿ_k⁻²⁾

(λ) (n

−

2)!

... ... ...

m k (λ) m

⁰

_k (λ) m k (λ)



 

  .

We are now in a position to define f(A) for arbitrary f!

For this let J A = diag(J 1 , J 2 , . . . , J k ) denote the Jordan canonical form of A, A = T J A T

⁻

¹ . Each Jordan block J j = J j (λ j ) is an (n j × n j ) matrix.

The characteristic polynomial c A of A is then given by c A (λ) = det(λI − A) =

Y k j=1

(λ − λ j ) ⁿ

^j

.

(9)

We now say that the function f is defined for A if

f is defined on an open set M containing the spectrum of A and f possesses n j − 1 derivatives at λ j .

In this case we define for j = 1, 2, . . . , k

f(J j (λ j )) :=



 



f(λ j ) f

⁰

(λ j ) · · · ^f

⁽^nj

_(n

⁻²⁾_j_−2)!

^(λ

^j

⁾ ^f

⁽^nj

_(n

⁻¹⁾_j_−1)!

^(λ

^j

⁾ f(λ j ) · · · ^f

(nj−3)

(λ

j

) (n

j−3)!

f

⁽^nj−²⁾

(λ

j

) (n

j−2)!

... ... ...

f (λ j ) f

⁰

(λ j ) f (λ j )



 



∈ C ⁿ

^j^×n^j

and, finally,

f(A) := T diag(f(J 1 ), f(J 2 ), . . . , f(J k )) T

⁻¹

.

30 Remarks

(1) For polynomials f(λ) = α 0 + α 1 λ + · · · + α m λ ^m the matrix f (A) thus defined is f(A) = α 0 I + α 1 A + · · · + α m A ^m , i.e. it coincides with the standard definition.

(2) An equivalent definition of f(A) uses the relation f(A) = p(A) for the Hermite interpolating polynomial to the values of f and it’s derivatives on the spectrum of A.

(3) For a rational function f = p/q with p ∈ P m and q ∈ P k , f is defined for A if, and only if, no pole of f coincides with an eigenvalue of A. In this case

f(A) = p(A)[q(A)]

⁻¹

= [q(A)]

⁻¹

p(A).

(4) All statements of Lemma 1.3 hold for any function f defined for A.

31 (5) If f is analytic in a neighborhood of 0 and has the Taylor

^∗

series f (λ) = P

_∞

j=0 α j λ ^j with convergence radius τ > 0 (τ = ∞ is allowed).

Then f is defined for any matrix A with spectral radius ρ(A) < τ and

f(A) = X

∞

j=0

α j A ^j = lim

m→∞

X m j=0

α j A ^j .

A well-known example for the last remark is the Neumann

^∗

Series (I − A)

⁻

¹ =

X

∞

j=0

A ^j , for ρ(A) < 1.

Brook Taylor (1685–1731)

(10)

The exponential function of a matrix A can, e.g., also be defined by exp(A) =

X

∞

j=0

1 j! A ^j . (6)

(The power series converges for any matrix A, since the (scalar) series has an infinite radius of convergence.)

The numerical determination of e ^A for A = diag( − 1, − 100) via (6) is prohibitive. Why?

Determine e ^A for the symmetric matrix A with eigenvalues − 1 and

− 100 and corresponding eigenvectors [1, 1] ^T and [1, − 1] ^T . Double precision is sufficient.

33 Theorem 1.4 (Properties of the matrix exponential).

Let A, B ∈ C ⁿ

^×

ⁿ , s, t ∈ R. Then

e ^(t+s)A = e ^tA e ^sA , e ^0A = I n (semigroup property), if AB = BA, then e ^A+B = e ^A e ^B ,

e ^tA is invertible with inverse (e ^tA )

⁻

¹ = e

⁻

^tA , t 7→ e ^tA is differentiable with _dt ^d e ^tA = Ae ^tA = e ^tA A.

In the analysis of ODEs it is often important to know the behavior of exp(tA) as t → ∞. The crucial quantity here is the spectral abscissa α(A) of A defined as

α(A) := max { Re(λ) : λ ∈ Λ(A) } .

Carl Gottfried Neumann (1832–1925)

34 Theorem 1.5 (Asymptotic behavior of exp(tA)).

Let A ∈ C ^n×n .

(1) There holds lim t

→∞

exp(tA) = 0 if, and only if, α(A) < 0.

(2) If α(A) > 0 then exp(tA) grows unbounded as t → ∞.

If α(A) = 0 then exp(tA) remains bounded as t → ∞ (but need not, in general, converge).

(3) There holds k exp(tA)k ≥ exp(tα(A)) for any matrix norm.

If A is normal (i. e. AA ^H = A ^H A), then k exp(tA) k 2 = exp(tα(A)).

For normal matrices k exp(tA)k 2 is thus a strictly decreasing function of t ≥ 0, if α(A) < 0.

If A is non-normal one observes the usual non-normality effects

(analogous to the behavior of kA ^m k 2 , m → ∞, if ρ(A) < 1).

(11)

Example

t 7→ k e ^tA

ⁱ

k 2 for the normal matrix A 1 =

₋

_{1 0}

0

−2

and the nonnormal A 2 =

₋

_{1 5}

0

−1

.

In both cases we have α(A i ) = − 1, but only for A 1 this function is exponentially decreasing, i. e. k e ^tA

ⁱ

k 2 = e

⁻

^t .

36 1.4 First Order Systems of ODEs

Only very few systems of the form (DE) or IVPs of the form (IVP) can be solved explicitly. Even first-order linear systems,

y

⁰

₁ = a 1,1 (t)y 1 + a 1,2 (t)y 2 + · · · + a 1,n (t)y n + b 1 (t), y

⁰

₂ = a 2,1 (t)y 1 + a 2,2 (t)y 2 + · · · + a 2,n (t)y n + b 2 (t),

... = ...

y

⁰

_n = a n,1 (t)y 1 + a n,2 (t)y 2 + · · · + a n,n (t)y n + b n (t) or, more briefly,

y

⁰

= A(t)y + b(t) with A(t) = [a i,j (t)] and b(t) = [b j (t)], (Lin) are among these exceptions only under further restrictions.

37 If the functions a i,j (t), b j (t) are continuous on an interval I and if k A(t) k ≤ L for all t ∈ I

(which we assume henceforth), then according to Theorem 1.1 (Lin) possesses a unique solution for any choice of initial condition

y

⁰

(t 0 ) = y 0 (t 0 ∈ I).

Satz 1.6 (Solutions of first-order linear systems).

The solutions of the homogeneous system y

⁰

= A(t)y form an n-dimensional subspace of C ¹ (I).

The difference of two solutions of the inhomogeneous system

y

⁰

= A(t)y + b(t) solves the associated homogeneous system.

(12)

In the special case of constant coefficients a i,j (t) ≡ a i,j for all t these solutions may be given in closed form. We consider first the

homogeneous case b(t) ≡ 0:

Let e 1 , e 2 , . . . , e n denote the unit coordinate vectors in R ⁿ . The functions

x j (t) := e ^(t−t

⁰

^)A e j , j = 1, 2, . . . , n, solve the IVP x

⁰

_j = Ax j , x j (t 0 ) = e j .

Moreover, the functions x 1 (t), x 2 (t), . . . , x n (t) are linearly independent and form a basis of the solution space of y

⁰

= Ay, i. e. each solution is of the form

y(t) = X n j=1

c j x j (t).

In particular, the solution for the IVP y

⁰

= Ay, y(t 0 ) = y 0 is given by y(t) = e ^(t

⁻

^t

⁰

^)A y 0 .

39 To solve the inhomogeneous IVP y

⁰

= Ay + b(t), y(t 0 ) = y 0 (for simplicity we assume the components b j (t) to be continuous on all of R), one may employ a technique known as variation of constants.

The (unique) solution is y(t) =

X n j=1

Z t t

0

w j (τ )

w(τ) dτ + y 0,j

x j (t), with the Wronski

^∗

determinants

w(t) = det X(t), w j (t) = det X j (t), where

X (t) =

x 1 (t) x 2 (t) . . . x n (t) X j (t) =

x 1 (t) · · · x j

−

1 (t) b(t) x j+1 (t) · · · x n (t)

∗

Josef Hoëné de Wronski (1778–1853)

40 An alternative representation formula for the solution of the inhomogeneous system with constant coefficients (also known by the name variation of constants) is

y(t) = exp((t − t 0 )A)y 0 + Z t

t

0

exp((t − τ )A) b(τ) dτ. (7) Note:

This formula also holds when b is a function of y as well as t, i.e.

b = b(t, y(t)).

(13)

The formula (7) also allows the analysis of the linearized problem.

Linearizing the differential equation y

⁰

= f(t, y) at the point (t 0 , y 0 ), one obtains (multivariate Taylor expansion)

f(t, y) ≈ f(t 0 , y 0 )

| {z }

=:b

+ f t (t 0 , y 0 )

| {z }

=:a

(t − t 0 ) + f

y

(t 0 , y 0 )

| {z }

=:A

(y − y 0 )

resulting in the linearized IVP

y

⁰

(t) = A(y − y 0 ) + (t − t 0 )a + b, y(t 0 ) = y 0 . (8) Applying (7) we obtain the solution of (8) as

y(t) = y 0 − (t − t 0 )A

⁻

¹ a + (e ^(t

⁻

^t

⁰

^)A − I)(A

⁻

¹ b + A

⁻

² a).

Confirm this formula by direct calculation or by a check of (8).

42 1.5 Further Examples

1.5.1 Satellite in the gravitation field of the Earth and Moon Assume the three bodies move in a fixed plane and that the two larger bodies rotate with constant separation distance and constant angular velocity around their common center of mass.

Thus satellite has no effect on the trajectories of Earth and Moon.

In a co-rotating coordinate system (with respect to which Earth and Moon are at rest) with their common center of mass as its origin, the satellite’s trajectory (x(t), y(t)) is described by

x

⁰⁰

= x + 2y

⁰

− µ

⁰

x + µ

[(x + µ) ² + y ² ] ^3/2 − µ x − µ

⁰

[(x − µ

⁰

) ² + y ² ] ^3/2 , y

⁰⁰

= y − 2x

⁰

− µ

⁰

y

[(x + µ) ² + y ² ] ^3/2 − µ y

[(x − µ

⁰

) ² + y ² ] ^3/2 .

43 µ = 1/82.45: ratio of Moon mass to total mass of Earth and Moon;

µ

⁰

= 1 − µ.

Unit of length is Earth-Moon-distance; Moon located on positive, Earth on negative x-axis.

Unit of time: angular velocity of rotation is one, i.e., Moon rotates around Earth in 2π time units.

Known as a restricted 3-body problem (because third body has no influence on remaining two).

Initial conditions: at time t = 0 satellite at position

(x(0), y(0)) = (1.2, 0) with velocity (x

⁰

(0), y

⁰

(0)) = (0, −1.05).

(14)

To transform this system into a first-order system set y 1 = x, y 2 = y, y 3 = x

⁰

, y 4 = y

⁰

, resulting in

y ₁

⁰

= y 3 , y ₂

⁰

= y 4 ,

y ₃

⁰

= y 1 + 2y 4 − µ

⁰

y 1 + µ

[(y 1 + µ) ² + y ₂ ² ] ^3/2 − µ y 1 − µ

⁰

[(y 1 − µ

⁰

) ² + y ₂ ² ] ^3/2 , y ₄

⁰

= y 2 − 2y 3 − µ

⁰

y 2

[(y 1 + µ) ² + y ₂ ² ] ^3/2 − µ y 2

[(y 1 − µ

⁰

) ² + y ₂ ² ] ^3/2 , with initial conditions

y 1 (0) = 1.2, y 2 (0) = y 3 (0) = 0 and y 4 (0) = − 1.05.

Solution (x(t), y(t)) = (y 1 (t), y 2 (t)) is a closed curve with period T ≈ 6.19.

45 Trajectory of the satellite

−1.5 −1 −0.5 0 0.5 1 1.5

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8

x(t)

y(t) Erde Mond

t=0.86132 t=1.4611

t=2.1265 t=3.0961

t=4.0657

t=4.731

t=5.3308

46 Components of velocity

0 1 2 3 4 5 6

−8

−6

−4

−2 0 2 4 6 8

t Geschwindigkeit in x−Richtung Geschwindigkeit in y−Richtung

(15)

1.5.2 Kinetics of an autocatalytic chemical reaction Three chemical substances S, T and U , autocatalytic reaction

S −→ ^k

¹

T, T + U −→ ^k

²

S + U, 2T −→ ^k

³

T + U.

Amounts of substances are described by the first-order ODE y

⁰

₁ = − k 1 y 1 + k 2 y 2 y 3 ,

y

⁰

₂ = k 1 y 1 − k 2 y 2 y 3 − 2k 3 y ₂ ² + k 3 y ₂ ² = k 1 y 1 − k 2 y 2 y 3 − k 3 y ² ₂ , y

⁰

₃ = k 3 y ² ₂ .

Reaction rates k j measure of speed at which the three reactions proceed. Typically orders of magnitude difference, e.g.

k 1 = 0.04, k 2 = 10 ⁴ , k 3 = 3 · 10 ⁷ .

48 For initial conditions y 1 (0) = 1; y 2 (0) = y 3 (0) = 0 we obtain:

0 200 400 600 800 1000

0 0.2 0.4 0.6 0.8 1

t yj(t)

Substanz S Substanz T Substanz U

Note:

Due to y

⁰

₁ (t) + y ₂

⁰

(t) + y ₃

⁰

In 1925 Umberto d’Ancona determines the percentage of sharks of total amount of fish caught in the harbor of Trieste:

1 Ordinary Differential Equations: Examples, theoretical foundations and tools

1.1 Volterra’s principle - a motivating example