3.

Integration and Simulation

Hil Meijer, Twente, NL
Outline

1. Phase Portraits; Local and Global – Matrix Exponents

2. Quadrature; Numerical Integration –Trapezoidal, Simpson
3. Simulations; Euler-Forward/Backward – Stability

Mathematical Methods 2023-11-28 2 / 29

Overview

Fitzhugh-Nagumo revisited
Phase Portraits
Diagonalization
Planar Linear Systems
Real Eigenvalues
Generating the portrait
Transitions

Quadrature - Numerical Integration

Simulations – Time stepping

Mathematical Methods 2023-11-28 3 / 29

FHN simulations

dV

dt = V (V − a)(1 − V ) − W + I,
dW
dt = ϵ(V − bW ).
Different parameter settings yield either
Subthreshold oscillations or Action potentials
0.1 1

0.8
0
0.6

0.4
−0.1
V

V
0.2

0
−0.2

−0.2

−0.3 −0.4
0 200 400 600 800 1000 0 200 400 600 800 1000
Time (ms) Time (ms)

Mathematical Methods 2023-11-28 4 / 29

FHN Phase Plane

Linearization yields a local approximating vector field.

Mathematical Methods 2023-11-28 5 / 29

From y ′ = Ay to Matrix Exponentials

Take the linear initial value problem

y ′ = f (t, y ) = Ay with y (t0 ) = y0 .
Look at y (t0 + t) = y(t0 ) + ty ′ + 12 t 2 y ′′ + 16 t 3 y ′′′ + ....
We need y ′′ and y ′′′ : y ′ = Ay, y ′′ = A2 y and y ′′′ = A3 y
In fact, we can show y (n) = An y .
P 
∞ t i Ai
This suggests the solution y (t0 + t) = i=0 i! y0 .
Proof: With A0 = In , we have y(t0 ) = y0 and also
∞ ∞ i i+1
! !
X t i−1 Ai X t A
y′ = y= y = Ay
(i − 1)! i!
i=0 i=0
∞ i i
X tA
We define the matrix exponent exp(At) :=
i!
i=0
Now we need to compute and understand this new matrix.
Mathematical Methods 2023-11-28 6 / 29
Diagonalization

Determine etA for  

λ1 0
A=
0 λ2
So
t k λk1
P∞ !
e λ1 t
 
tA k=0 k! 0 0
e = =
0
P∞ t k λk2 0 e λ2 t
k=0 k!

Try a linear, invertible coordinate transformation y = Tu:

u ′ = T −1 y ′ = T −1 Ay = T −1
| {zAT} u
D
If we choose T to be the matrix with eigenvectors as columns,
then D is diagonal. (T maps standard unit vectors to eigenvectors)
Mathematical Methods 2023-11-28 7 / 29
Diagonalization

Determine etA for

   k 
λ1 0 k λ1 0
A= =⇒ A =
0 λ2 0 λk2
So
t k λk1
P∞ !
e λ1 t
 
tA k=0 k! 0 0
e = =
0
P∞ t k λk2 0 e λ2 t
k=0 k!

Try a linear, invertible coordinate transformation y = Tu:

u ′ = T −1 y ′ = T −1 Ay = T −1
| {zAT} u
D
If we choose T to be the matrix with eigenvectors as columns,
then D is diagonal. (T maps standard unit vectors to eigenvectors)
Mathematical Methods 2023-11-28 7 / 29
Eigenvalue combinations

 
′ a11 a12
y = y
a21 a22
Get eigenvalues from characteristic polynomial

p(λ) = det(A − λI) = λ2 − (a11 + a22 ) λ + (a11 a22 − a21 a12 )

| {z } | {z }
T D
2
= λ − Tλ + D = 0
Possible eigenvalue combinations (Consider sign T 2 − 4D)
▶ two different real eigenvalues T 2 − 4D > 0
▶ a complex pair T 2 − 4D < 0
▶ two real, equal eigenvalues T 2 − 4D = 0 (ignored in this
course)
Mathematical Methods 2023-11-28 8 / 29
Real Eigenvalues

Possibilities
▶ λ1 < 0 < λ2 Saddle
▶ λ1 < λ2 < 0 Stable Node
▶ 0 < λ1 < λ2 Unstable Node (reverse time)
General solution is x(t) = c1 eλ1 t v1 + c2 eλ2 t v2

For c1 = 0 or c2 = 0 the semi-axes are solutions (orbits) of the

system. These are invariant.
Mathematical Methods 2023-11-28 9 / 29
Phase Portrait for Stable Node λ1 < λ2 < 0

Solution form: x(t) = c1 eλ1 t v1 + c2 eλ2 t v2

▶ If t → +∞, then solutions become parallel to v2 .
▶ If t → −∞, then solutions become parallel to v1 .
▶ If v1 = (1, 0)T en v2 = (0, 1)T , then xi = ci eλi t .
 λ2 /λ1
log x1 /c1 log x2 /c2 x1
=t = =⇒ x2 = c2
λ1 λ2 c1
Orbits in (v1 , v2 )-basis look like k(= λ2 /λ1 )-power
functions.

Mathematical Methods 2023-11-28 10 / 29

Summary Phase Portraits for Real Eigenvalues
Recipe for drawing:
▶ Sketch the lines spanned by the eigenvectors.
▶ Add arrows according to real part of eigenvalue: If λ < 0
then towards origin
▶ Add orbits just off the lines spanned by the eigenvectors
▶ Saddle: Start along the stable direction and move away
from the origin
▶ Node: Assuming |λ2 | < |λ1 |. Orbits are parallel to
eigenvector v2 near origin and parallel to v1 if far away.

Mathematical Methods 2023-11-28 11 / 29

Complex eigenvalues: Phase Portrait

Recipe for drawing:

▶ Eigenvectors are of no use.
▶ When ℜ(λ) = α < 0, orbits spiral towards
the origin.
▶ Spirals or ellips are horizontal when
x2′ = 0 and vertical when x1′ = 0.
▶ Direction of Rotation: Compute the vector
field at a convenient point, e.g. (1, 0). Use
that to see whether it is clockwise or
counterclockwise.

Mathematical Methods 2023-11-28 12 / 29

Generating a phase portrait

%Time array and system y’=Ay

t=-1:.05:7;alpha=2; %Vary -8<alpha<+6
A=[-2 4; alpha -5];
%Initial conditions
y0=[1 1 1 0 0 -1 -1 -1; 1 0 -1 1 -1 1 0 -1];
figure(1);clf(1);hold on;
plot(y0(1,:),y0(2,:),’LineStyle’,none,’Marker’,’o’);
for i=1:8;
for jj=1:size(t,2)
sol(1:2,jj)=expm(t(jj)*A)*y0(:,i);
end
plot(sol(1,:),sol(2,:));
end
xlim([-2 2]);ylim([-2 2]);

Mathematical Methods 2023-11-28 13 / 29

Improving the code

▶ Axes should have labels. Always.

▶ wrong coordinates selection to plot, use
plot(y0(:,1),y0(:,2)
▶ Pre-allocation of memory of solution vector.
▶ The order of the for-loops is wrong; expm is computed for
every orbit but rather you compute it once for every t and
then use it for every initial condition. Changing dt=0.05 to
dt=0.001, you will feel the difference.

Mathematical Methods 2023-11-28 14 / 29

Generating a phase portrait – Improved

%Time array and system y’=Ay

t=-1:.05:7;alpha=2; %Vary -8<alpha<+6
A=[-2 4; alpha -5];
%Initial conditions
y0=[1 1 1 0 0 -1 -1 -1; 1 0 -1 1 -1 1 0 -1];
figure(1);clf(1);hold on;
plot(y0(:,1),y0(:,2),’LineStyle’,none,’Marker’,’o’);
sol=nan(2,length(t),8);
for jj=1:length(t)
EA=expm(t(jj)*A);
for i=1:8;
sol(1:2,jj,i)=EA*y0(:,i);
end
end
for i=1:8
plot(sol(1,:,i),sol(2,:,i),’Linewidht’,1);
end
xlim([-2 2]);ylim([-2 2]);
xlabel(’X_{1}’);ylabel(’X_{2}’,’Rotation’,0);

Mathematical Methods 2023-11-28 15 / 29

Translating to Python with ChatGPT

https://chat.openai.com/share/
07e3ac2a-9e2d-4347-a948-5b8160c8bb4b
The code appears to be pretty ok. Execution leads to two
problems.
▶ First, The command Hold is not there, and is not
needed=⇒ Remove that line
▶ Second, The command expm is part of the package
scipy.linalg not numpy. =⇒ Load that submodule explicitly.

Mathematical Methods 2023-11-28 16 / 29

Translating to Python with ChatGPT
import numpy as np
import scipy.linalg as spla
import matplotlib.pyplot as plt
# Time array and system y’=Ay
t = np.arange(-1, 7.05, 0.05)
alpha = 2 # Vary -8 < alpha < +6
A = np.array([[-2, 4], [alpha, -5]])
# Initial conditions
y0 = np.array([[1, 1, 1, 0, 0, -1, -1, -1],
[1, 0, -1, 1, -1, 1, 0, -1]])
plt.figure(1)
plt.plot(y0[0, :], y0[1, :], linestyle=’none’, marker=’o’)

for i in range(8):
sol = np.zeros((2, len(t)))
for jj in range(len(t)):
sol[:, jj] = np.dot(spla.expm(t[jj] * A), y0[:, i])
plt.plot(sol[0, :], sol[1, :])

plt.xlim([-2, 2])
plt.ylim([-2, 2])
plt.show()
Mathematical Methods 2023-11-28 17 / 29
The transition from saddle to spiral
saddle critical node stable spiral

α=2 α = 9/2 α=5

Mathematical Methods 2023-11-28 18 / 29

Overview

Fitzhugh-Nagumo revisited
Phase Portraits
Diagonalization
Planar Linear Systems
Real Eigenvalues
Generating the portrait
Transitions

Quadrature - Numerical Integration

Simulations – Time stepping

Mathematical Methods 2023-11-28 19 / 29

Quadrature

Quadrature:=Numerical integration
Z b
f (x)dx =?
a

Section 19.5 on Trapezoidal rule and Simpson rule

Mathematical Methods 2023-11-28 20 / 29

Trapezoidal Rule: The Algorithm

Divide the interval [a,b] in N subintervals of length h so

b − a = Nh and xi = a + h ∗ i.
Approximate the local area by averaging the function using the
endpoints of the subinterval:
First interval: Area = width × height = h × 12 (f (a) + f (a + h))
Second interval: Area = h × 21 (f (a + h) + f (a + 2h))
Summing up:
Rb
(f (a) + f (b)) + h N
h P
a f (x)dx ≈ 2 i=1 f (a + ih)
 
1 1
=h 1 ... 1 .f (⃗x )
2 2
| {z }
w

Mathematical Methods 2023-11-28 21 / 29

Trapezoidal Rule: error estimate

′′
For one interval: f (x) − p1 (x) = (x − x0 )(x − x1 ) f2
Determine error by integration
Rx Rx ′′
E1 = x01 f (x) − p1 (x)dx = x01 (x − x0 )(x − x1 ) f2 dx
f ′′ h3 ′′
Rh
|{z} 0 v (v − h) 2 dv = 12 f
=
x=x0 +v

Take C = maxx∈[a,b] (f ′′ (x)) and sum over all intervals:

(b − a)h2
Error ≈ C
12

Mathematical Methods 2023-11-28 22 / 29

Simpson Rule: Algorithm

Divide the interval [a,b] in 2N subintervals of length h so b − a = Nh

Fit a quadratic polynomial through 3 subsequent points:

(x − x2 )(x − x1 ) (x − x2 )(x − x0 ) (x − x0 )(x − x1 )

p2 (x) = f (x0 ) + f (x1 ) + f (x2 )
(x0 − x2 )(x0 − x1 ) (x1 − x2 )(x1 − x0 ) (x2 − x0 )(x2 − x1 )

Show that xx2 p2 (x)dx = h

R
3
(f0 + 4f1 + f2 )
0
Summing up:
Rb h
a f (x)dx ≈ 3
(f0 + 4f1 + 2f2 + 4f3 + · · · + 2f2N−2 + 4f2N−1 + f2N )
h
= (1 4 2 4 2 ... 2 4 1) .f (⃗x )
3
| {z }
w

With fi = f (x + ih) and error estimate

ε = Ch4

for some constant C.

Mathematical Methods 2023-11-28 23 / 29

Overview

Fitzhugh-Nagumo revisited
Phase Portraits
Diagonalization
Planar Linear Systems
Real Eigenvalues
Generating the portrait
Transitions

Quadrature - Numerical Integration

Simulations – Time stepping

Mathematical Methods 2023-11-28 24 / 29

Approximating ODEs

Nonlinear, nonautonomous ODE y ′ = f (x, y ):

▶ No explicit solution, but approximate it numerically at
discrete time-steps.
▶ Algorithm can only use function evaluations.
▶ Equal timesteps: non-adaptive
▶ Changing timesteps: adaptive, depending on local error
estimates.

Mathematical Methods 2023-11-28 25 / 29

Euler Forward

Start (again) with y ′ = f (x, y ), but now do not assume

fx (x, y ) = 0.
(1) Look for a solution y(x) (2) Take a step h in time (3) Taylor
again
y(x + h) = y(x) + hy ′ + O(h2 ) = y(x) + hf (x, y) + O(h2 )
Ignore the quadratic terms of h. So, we get an numerical
approximation of the solution at the next step:
y (x) = y0 , y1 = y(x+h) = y0 +hf (x, y0 ), y(x+2h) = y1 +hf (x+
Error estimates O(h2 ) = 21 h2 y ′′ (x̃)
We take (xend − xstart )/h steps and at each the error is
proportional to h2 :
global error: proportional to h.
Mathematical Methods 2023-11-28 26 / 29
Euler Backward

Integrating dy = f (x, y )dx from x to x + h yields

Z x+h
y(x + h) − y(x) = f (s, y(s))ds.
x

The approximation
▶ f (s, y (s)) ≈ f (x, y(x)) yields Euler-forward.
▶ f (s, y (s)) ≈ f (x + h, y (x + h)) yields Euler-backward.
Euler-backward defines the solution implicitly by

y (x + h) − y (x) = hf (x + h, y (x + h)).

Finding a zero, e.g. Newton-iteration, at each time-step

needed.

Mathematical Methods 2023-11-28 27 / 29

Euler Forward: Stability

Consider the simple testproblem y ′ = λy.

Start at y (0) = y0 and take steps of length h. Euler Forward
generates iteration given by
yn = yn−1 + hf (xn−1 , yn−1 ) = yn−1 + hλyn−1 = (1 + hλ)yn−1
If Re(λ) < 0, then y (t) → 0 as t → ∞.
Stability of yn = (1 + hλ)n y0 implies |1 + hλ| < 1.
Note that 1 + hλ is the eigenvalue of the fixed point 0 of the
iteration.

Mathematical Methods 2023-11-28 28 / 29

Euler Backward: Stability

Consider the simple testproblem y ′ = λy.

Start at y (0) = y0 and takes steps of length h. Euler Backward
generates iteration implicitly given by
yn = yn−1 + hf (xn , yn ) = yn−1 + hλyn
Rearranging terms yields
1
yn = yn−1
1 − hλ
If Re(λ) < 0, then y (t) → 0 as t → ∞.
Stability of yn = (1 − hλ)−n y0 implies |1 − hλ| > 1.

Mathematical Methods 2023-11-28 29 / 29