Lecture 7 · ODE and Simulation

Lecture 7 · ODE and Simulation · October 21

ODE and Simulation

How does a simulator advance the world?

$$\dot x=f(x,t),\qquad x(t_0)=x_0$$

x: state column
f: rate of change
t: time; t₀ is the start
x₀: known initial state

Ordinary differential equation (ODE): a state’s rate of change with time.

Mathematics for Robotics · ROB-GY 6013 · NYU Tandon · Jing Zhang · Fall 2026

Initial-value problem (IVP) — [4] FNC §6.1

Objects

ODE and Simulation

A trajectory and its stored samples

$$t_k=t_0+kh,\qquad x_k\approx x(t_k)$$

h, k: step size and sample index
x(t), xₖ: exact state and numerical sample

Initial state + rate model

↓ exact solution / numerical integration

Trajectory at requested times

Accuracy: error · Stability: repeated steps

Solver verification: exact / refined solution

Model validation: measurements

Exact solution; numerical solution — [4] FNC §6.2 · Course map

Definition

Vector Fields

An arrow gives a rate, not a destination

$$\dot x=f(x,t)\in\mathbb R^d$$

x: d-component state column
f(x,t): d-component state rate
t: time held fixed in a field picture

Illustration: x=(q,v)ᵀ and f=(v,−sin q)ᵀ. Arrow lengths are scaled.

Vector field; phase portrait — [1] Tedrake, ch.2 · Course illustration

Assumptions

Vector Fields

One field, different initial states

$$\dot x=-2x,\quad x(0)=1\ \text{or}\ 2$$

x(0): two possible starts
−2x: the same rate model for both

Initial-value problem (IVP): a differential equation plus a starting state.

Initial condition; solution trajectory — [4] FNC §6.1 · Course example

Worked example

Vector Fields

Verify an exact decay solution

$$x(t)=e^{-2t}\quad\Longrightarrow\quad\dot x=-2x$$

x(0): one
−2: decay rate in inverse time
t: elapsed time

Worked example: x(0.1) = 0.8187307531.

Exponential decay — [4] FNC §6.1 · Course example

Your turn

Vector Fields

Your turn: a changed decay rate

$$\dot x=-3x,\qquad x(0)=2$$

x: scalar state
−3: fixed decay rate
2: initial value

Which trajectory satisfies BOTH the rate and the start?

exp(−3t)
2 exp(−3t)
2 exp(3t)

✓

Correct B. Check start AND rate.

Solution 1 / 4

1 · Solve the rate equation

x(t): scalar state at time t · c: constant set by the initial value

For a constant decay rate −3, the exponential family is $$x(t)=ce^{-3t},\qquad \dot x(t)=-3ce^{-3t}=-3x(t).$$

Solution 2 / 4

2 · Impose the initial value

x(t): scalar state at time t · c: constant set by the initial value

At time zero the exponential is one, so the constant must be two. $$x(0)=ce^0=c=2\quad\Longrightarrow\quad x(t)=2e^{-3t}.$$

Solution 3 / 4

3 · Verify both conditions

x(t): scalar state at time t · c: constant set by the initial value

The chosen trajectory starts at two and has the required derivative. $$x(0)=2,\qquad \dot x(t)=-6e^{-3t}=-3x(t).$$

Solution 4 / 4

Answer · B

x(t): scalar state at time t · c: constant set by the initial value

A starts at one, not two.
B has the correct initial value and decay rate.
C reverses the sign of the rate.

Exponential decay — [4] FNC §6.1 · Course example

Why

Robot State

Position alone does not predict motion

$$x=\begin{bmatrix}q\\v\end{bmatrix},\qquad v=\dot q$$

q: n joint positions
v: n joint velocities
x: 2n-component state

The same angle can occur while a joint moves in either direction.

First-order system; state vector — [4] FNC §6.3

Model

Robot State

The manipulator equation is a force balance

$$M(q)\ddot q+C(q,v)v+g(q)=\tau$$

M: n×n inertia; zᵀMz>0 for every z≠0
Cv: velocity-dependent generalized forces
g, τ: gravity and applied generalized forces
q,v: positions and velocities

The model balances inertial, velocity-dependent and gravitational terms.

Manipulator dynamics; mass matrix — [2] Modern Robotics §8.1

Derivation

Robot State

Rewrite the mechanics as a first-order ODE

$$\dot x=\begin{bmatrix}v\\M(q)^{-1}[\tau-C(q,v)v-g(q)]\end{bmatrix}$$

x: [q;v], a 2n-column
τ: applied generalized force
M⁻¹: exists for positive-definite inertia

First order: only ẋ appears. Solve M a = τ − Cv − g for acceleration.

First-order system; forward dynamics — [2] Modern Robotics §8.1 · [4] FNC §6.3

Worked example

Robot State

A one-joint calculation

$$a=\frac{5-0.5(2)-3}{2}=0.5,\qquad\dot x=\begin{bmatrix}2\\0.5\end{bmatrix}$$

M,C: 2 and 0.5 at this state
g,τ: 3 and 5
q,v: 1 and 2

Worked example: solve acceleration, then assemble the full derivative.

Forward dynamics — [2] Modern Robotics §8.1 · Course example

Your turn

Robot State

Your turn: reduce the applied force

$$\begin{gathered}M=2,\ C=.5,\ g=3,\ \tau=1,\ v=2\\ Ma=\tau-Cv-g,\quad\dot x=(v,a)^T\end{gathered}$$

a: acceleration
x=(q,v)ᵀ: position then velocity

Which state rate follows from the force balance?

(−1.5,2)ᵀ
(2,−3)ᵀ
(2,−1.5)ᵀ

✓

Correct C. Velocity, then acceleration.

Solution 1 / 5

1 · Compute the velocity force

q,v,a: position, velocity and acceleration · M,C: inertia and velocity-force coefficient · τ,g: applied force and gravity force

The given coefficient is 0.5 and the velocity is 2. $$Cv=0.5\times2=1.$$

Solution 2 / 5

2 · Apply force balance

q,v,a: position, velocity and acceleration · M,C: inertia and velocity-force coefficient · τ,g: applied force and gravity force

Subtract the velocity force and gravity from the applied force. $$Ma=\tau-Cv-g=1-1-3=-3.$$

Solution 3 / 5

3 · Divide by inertia

q,v,a: position, velocity and acceleration · M,C: inertia and velocity-force coefficient · τ,g: applied force and gravity force

Acceleration is net force divided by inertia; it is not the net force itself. $$a=\frac{-3}{2}=-1.5.$$

Solution 4 / 5

4 · Keep the state order

q,v,a: position, velocity and acceleration · M,C: inertia and velocity-force coefficient · τ,g: applied force and gravity force

Position differentiates to velocity; velocity differentiates to acceleration. $$x=(q,v)^T\quad\Longrightarrow\quad\dot x=(v,a)^T=(2,-1.5)^T.$$

Solution 5 / 5

Answer · C

q,v,a: position, velocity and acceleration · M,C: inertia and velocity-force coefficient · τ,g: applied force and gravity force

A swaps velocity and acceleration.
B omits division by inertia two.
C preserves component order and divides net force by inertia.

Forward dynamics — [2] Modern Robotics §8.1 · Course example

Model

Linear Dynamics

A constant matrix gives a linear field

$$\dot x=Ax,\qquad A\in\mathbb R^{d\times d}$$

A: constant matrix
x: d-component state
Ax: state rate at x

Illustrated A has rows (−1,−1) and (0,−2); the field is homogeneous.

Linear constant-coefficient system — [3] Strang, lecture 23 · Course example

Recall

Linear Dynamics

Recall the eigenvector coordinates

$$AV=V\Lambda,\qquad x=Vc$$

V: d independent eigenvectors as columns
Λ: eigenvalues in matching order
c: coordinates of x in that basis

Eigenbasis: d independent eigenvectors spanning the state space. Here v₁=(1,0)ᵀ, v₂=(1,1)ᵀ.

Diagonalization; modal coordinates — [3] Strang, lecture 23

Derivation

Linear Dynamics

The ODE separates into scalar modes

$$V\dot c=AVc=V\Lambda c\quad\Rightarrow\quad\dot c=\Lambda c$$

V: constant invertible basis
cᵢ: coefficient of eigenvector i; one mode
Λ: diagonal rate matrix

Multiply on the left by V⁻¹; each coordinate now evolves alone.

Decoupled linear ODEs — [3] Strang, lecture 23 · Course derivation

Result

Linear Dynamics

Each eigenvalue is a continuous growth rate

$$c_i(t)=e^{\lambda_i t}c_i(0)$$

λᵢ: rate of mode i
cᵢ(0): initial amount of that mode
t: time measured from zero

Negative rates decay; positive rates grow if the mode is present.

Eigenvalues; exponential growth and decay — [3] Strang, lecture 23 · Course example

Your turn

Linear Dynamics

Return from modal coordinates

$$V=\begin{bmatrix}1&1\\0&1\end{bmatrix},\quad \Lambda=\operatorname{diag}(-1,-2)$$

V: eigenvector basis
Λ: matching eigenvalues
c₀: coordinates of x₀ in the V basis

For NEW x₀=(1,1)ᵀ, which modal description is correct?

c₀=(0,1)ᵀ; only the −2 mode
c₀=(1,2)ᵀ; both modes
c₀=(1,1)ᵀ; both modes

✓

Correct A. Use the inverse basis.

Solution 1 / 5

1 · Invert the coordinate map

x₀: new initial state (1,1)ᵀ · V,Λ: eigenvector matrix and matching eigenvalues · c₀,c(t): initial and evolving modal coordinates

Because x=Vc, solve for c using the inverse, not the transpose. $$V=\begin{bmatrix}1&1\\0&1\end{bmatrix},\quad V^{-1}=\begin{bmatrix}1&-1\\0&1\end{bmatrix}.$$

Solution 2 / 5

2 · Transform the new start

x₀: new initial state (1,1)ᵀ · V,Λ: eigenvector matrix and matching eigenvalues · c₀,c(t): initial and evolving modal coordinates

Multiply each row of the inverse by the initial column. $$c_0=V^{-1}\begin{bmatrix}1\\1\end{bmatrix}=\begin{bmatrix}1-1\\1\end{bmatrix}=\begin{bmatrix}0\\1\end{bmatrix}.$$

Solution 3 / 5

3 · Evolve the active mode

x₀: new initial state (1,1)ᵀ · V,Λ: eigenvector matrix and matching eigenvalues · c₀,c(t): initial and evolving modal coordinates

The coefficient of the −1 mode is zero; only the −2 mode remains. $$c(t)=\begin{bmatrix}0e^{-t}\\1e^{-2t}\end{bmatrix}=\begin{bmatrix}0\\e^{-2t}\end{bmatrix}.$$

Solution 4 / 5

4 · Return to state coordinates

x₀: new initial state (1,1)ᵀ · V,Λ: eigenvector matrix and matching eigenvalues · c₀,c(t): initial and evolving modal coordinates

Both state components are equal even though only one modal coordinate is active. $$x(t)=Vc(t)=\begin{bmatrix}e^{-2t}\\e^{-2t}\end{bmatrix},\qquad x(0)=\begin{bmatrix}1\\1\end{bmatrix}.$$

Solution 5 / 5

Answer · A

x₀: new initial state (1,1)ᵀ · V,Λ: eigenvector matrix and matching eigenvalues · c₀,c(t): initial and evolving modal coordinates

A recognizes the second eigenvector.
B uses the transpose instead of the inverse of a nonorthogonal basis.
C copies state coordinates without changing basis.

Change of basis; matrix exponential — [3] Strang, lecture 23 · Course derivation

Definition

Linear Dynamics

The matrix exponential always exists

$$e^{At}=\sum_{j=0}^{\infty}\frac{(At)^j}{j!},\qquad x(t)=e^{At}x_0$$

A: any square constant matrix
j: series index; A⁰=I
x₀: initial state

Diagonalization is one way to evaluate the solution, not an existence requirement.

Matrix exponential — [3] Strang, lecture 23

Your turn

Linear Dynamics

A defective system

$$A=\begin{bmatrix}-1&1\\0&-1\end{bmatrix},\quad x_0=(0,1)^T,\quad x(t)=(te^{-t},e^{-t})^T$$

defective: fewer than n independent eigenvectors in ℂⁿ
v: eigenvector; here v lies along (1,0)ᵀ

Does x(t) satisfy ẋ=Ax and x(0)=x₀?

No: a defective matrix has no solution.
No: the factor t forces growth forever.
Yes: its derivative and initial state both match.

✓

Correct C. Both conditions hold.

Solution 1 / 5

1 · Differentiate the candidate

A: given constant state matrix · x(t),x₀: candidate trajectory and required initial state · t: continuous time

Use the product rule for the first component. $$\frac{d}{dt}(te^{-t})=e^{-t}-te^{-t}=(1-t)e^{-t},\quad\frac{d}{dt}e^{-t}=-e^{-t}.$$

Solution 2 / 5

2 · Compute the model rate

A: given constant state matrix · x(t),x₀: candidate trajectory and required initial state · t: continuous time

Multiply A by the candidate and compare component by component. $$Ax=\begin{bmatrix}-1&1\\0&-1\end{bmatrix}\begin{bmatrix}te^{-t}\\e^{-t}\end{bmatrix}=\begin{bmatrix}(1-t)e^{-t}\\-e^{-t}\end{bmatrix}=\dot x.$$

Solution 3 / 5

3 · Check the start

A: given constant state matrix · x(t),x₀: candidate trajectory and required initial state · t: continuous time

The candidate satisfies the initial value as well as the differential equation. $$x(0)=\begin{bmatrix}0e^0\\e^0\end{bmatrix}=\begin{bmatrix}0\\1\end{bmatrix}=x_0.$$

Solution 4 / 5

4 · Interpret the solution

A: given constant state matrix · x(t),x₀: candidate trajectory and required initial state · t: continuous time · v: eigenvector · I: identity matrix

Here λ=−1 is repeated. Its eigenvectors have only one independent direction. $$(A+I)v=0\ \Rightarrow\ v_2=0\ \Rightarrow\ v=c(1,0)^T,\quad c\ne0.$$
The exponential eventually dominates t. The matrix exponential gives the solution even without an eigenbasis. $$\lim_{t\to\infty}te^{-t}=0,\qquad x(t)=e^{At}x_0.$$

Solution 5 / 5

Answer · C

A: given constant state matrix · x(t),x₀: candidate trajectory and required initial state · t: continuous time

A confuses an incomplete eigenbasis with absence of a solution.
B ignores the exponential decay multiplying t.
C checks both defining conditions.

Defective matrix; matrix exponential — [3] Strang, lecture 23 · [7] NIST DLMF §1.2 · Course example

Worked example

Linear Dynamics

A real rotation has complex eigenvalues

$$A=\begin{bmatrix}0&-1\\1&0\end{bmatrix},\quad x(t)=(\cos t,\sin t)^T$$

x₀: (1,0)ᵀ
A: constant rotation generator
λ: ±i over the complex field

This trajectory stays on the unit circle.

Complex eigenvalues; rotation — [3] Strang, lecture 23 · Course example

Your turn

Linear Dynamics

Zero real parts need a closer look

$$A=\begin{bmatrix}0&1\\0&0\end{bmatrix},\qquad x_0=(0,1)^T$$

A: both eigenvalues zero
x₀: chosen initial state
x(t): state to determine

What does this zero-eigenvalue example do?

x(t)=(0,1)ᵀ: fixed
x(t)=(t,1)ᵀ: unbounded
x(t)=(exp(t),1)ᵀ: exponential growth

✓

Correct B. Coupling gives linear growth.

Solution 1 / 5

1 · Read the component equations

x₁(t),x₂(t): first and second coordinates; subscripts here label coordinates · A: given matrix coupling coordinate 2 into coordinate 1 · s: integration variable

Matrix multiplication gives one constant coordinate and one driven coordinate. $$\dot x=Ax\quad\Longrightarrow\quad\dot x_1=x_2,\quad\dot x_2=0.$$

Solution 2 / 5

2 · Solve the second coordinate

x₁(t),x₂(t): first and second coordinates; subscripts here label coordinates · A: given matrix coupling coordinate 2 into coordinate 1 · s: integration variable

Zero derivative means constant; its initial value is one. $$x_2(t)=x_2(0)+\int_0^t0\,ds=1.$$

Solution 3 / 5

3 · Integrate the first coordinate

x₁(t),x₂(t): first and second coordinates; subscripts here label coordinates · A: given matrix coupling coordinate 2 into coordinate 1 · s: integration variable

The first coordinate starts at zero and has constant derivative one. $$x_1(t)=x_1(0)+\int_0^t1\,ds=0+t=t.$$

Solution 4 / 5

4 · Classify the result

x₁(t),x₂(t): first and second coordinates; subscripts here label coordinates · A: given matrix coupling coordinate 2 into coordinate 1 · s: integration variable

The state grows linearly, despite both eigenvalues being zero. $$x(t)=(t,1)^T,\qquad\|x(t)\|_2=\sqrt{t^2+1}\to\infty.$$

Solution 5 / 5

Answer · B

x₁(t),x₂(t): first and second coordinates; subscripts here label coordinates · A: given matrix coupling coordinate 2 into coordinate 1 · s: integration variable

A ignores coupling between coordinates.
B integrates the constant second coordinate into linear growth.
C invents an exponential where growth is linear.

Boundedness; asymptotic decay — [3] Strang, lecture 23

Derivation

Euler's Method

An exact step is an integral

$$x(t+h)-x(t)=\int_t^{t+h}f(x(s),s)\,ds$$

h: positive interval length
s: integration variable
f(x(s),s): rate along the unknown path

The fundamental theorem of calculus gives the starting relation.

Integral form of an ODE — [4] FNC §6.2 · Course derivation

Rule

Euler's Method

Euler uses the rate at the start

$$x_{k+1}=x_k+h f(x_k,t_k)$$

xₖ: current numerical state
h: positive time step
f(xₖ,tₖ): rate at the old state

Every component uses the same old state.

Explicit Euler method — [4] FNC §6.2

Worked example

Euler's Method

One step against an exact answer

$$x_1=1+0.1(-2)=0.8$$

x₀: one
h: 0.1
f(x): −2x

Worked example: exact x(0.1)=0.8187307531; Euler is lower.

Explicit Euler method; one-step error — [4] FNC §6.2 · Course example

Your turn

Euler's Method

Your turn: two larger steps

$$\dot x=-2x,\quad x_0=1,\quad h=0.2$$

h: fixed time step
x₀: initial numerical state
x₂: approximation at t=0.4

Which Euler pair uses a fresh rate at each step?

x₁=.6, x₂=.36
x₁=.6, x₂=.2
x₁=.8, x₂=.64

✓

Correct A. Reevaluate the rate.

Solution 1 / 4

1 · Write the Euler update

xₖ: accepted numerical state after k steps · h,tₖ: step size 0.2 and time tₖ=kh · f(x): state rate −2x

Each step uses the rate at the current accepted state. $$x_{k+1}=x_k+hf(x_k)=x_k+0.2(-2x_k)=0.6x_k.$$

Solution 2 / 4

2 · Take the first step

xₖ: accepted numerical state after k steps · h,tₖ: step size 0.2 and time tₖ=kh · f(x): state rate −2x

At x₀=1, the rate is −2. $$f(x_0)=-2,\qquad x_1=1+0.2(-2)=0.6.$$

Solution 3 / 4

3 · Recompute the second rate

xₖ: accepted numerical state after k steps · h,tₖ: step size 0.2 and time tₖ=kh · f(x): state rate −2x

The current state is now 0.6, so the rate changes to −1.2. $$f(x_1)=-2(0.6)=-1.2,\qquad x_2=0.6+0.2(-1.2)=0.36.$$

Solution 4 / 4

Answer · A

xₖ: accepted numerical state after k steps · h,tₖ: step size 0.2 and time tₖ=kh · f(x): state rate −2x

A reevaluates the rate at the new state.
B freezes the initial rate for both steps.
C uses half the given step size.

Explicit Euler method — [4] FNC §6.2 · Course example

Error analysis

Euler's Method

Local error and accumulated error differ

$$\begin{aligned}x(t+h)&=x+h\dot x+\tfrac12h^2\ddot x+O(h^3)\\\ddot x&=f_t+J_f f\end{aligned}$$

x, f: evaluated on the exact path at time t
fₜ, Jf: time derivative; (Jf)ᵢⱼ=∂fᵢ/∂xⱼ

C²: continuous derivatives through order 2. O(hᵖ): size ≤ C hᵖ for small h. Local state error = h × normalized local truncation error.

Local error; global error — [4] FNC §6.2

Recall and derivation

Numerical Stability

A repeated step has a multiplier

$$\dot x=\lambda x\quad\Rightarrow\quad x_{k+1}=(1+h\lambda)x_k$$

λ: constant scalar rate
h: step size
1+hλ: numerical amplification factor
H: Hessian: second derivatives of the fixed quadratic

Recall Lecture 6: A=−H turns this into the gradient-descent multiplier.

Amplification factor — [4] FNC §11.3

Derivation

Numerical Stability

The strict inequality decides decay

$$|1+h\lambda|<1\quad\Longleftrightarrow\quad0<h<-2/\lambda$$

λ: real negative rate
h: positive step size
|1+hλ|: size of each numerical multiplier

FNC: bounded if |1+hλ| ≤ 1. Here: decay needs |1+hλ| < 1.

Absolute stability; asymptotic decay — [4] FNC §11.3

Your turn

Numerical Stability

Predict three step-size regimes

$$\dot x=-4x,\quad h\in\{.25,.4,.6\}$$

λ: fixed decay rate
h: candidate steps
x₀: any nonzero start

For h=.25, .4, .6, which Euler behaviors occur in order?

Zero; alternating growth; alternating growth
All decay because the ODE is stable
Zero; alternating decay; alternating growth

✓

Correct C. Magnitude decides decay.

Solution 1 / 5

1 · Derive the multiplier

h: positive time step · a: Euler multiplier · xₖ: accepted numerical state; x₀≠0

Its sign controls sign changes; its magnitude controls decay or growth. $$x_{k+1}=x_k+h(-4x_k)=ax_k,\quad a=1-4h,\quad x_k=a^kx_0.$$

Solution 2 / 5

2 · Test h = 0.25

h: positive time step · a: Euler multiplier · xₖ: accepted numerical state; x₀≠0

A zero multiplier sends the state to zero in one step. $$a=1-4(0.25)=0\quad\Longrightarrow\quad x_1=0.$$

Solution 3 / 5

3 · Test h = 0.4

h: positive time step · a: Euler multiplier · xₖ: accepted numerical state; x₀≠0

The negative sign alternates; the magnitude below one makes the states decay. $$a=1-4(0.4)=-0.6,\qquad |a|=0.6<1.$$

Solution 4 / 5

4 · Test h = 0.6

h: positive time step · a: Euler multiplier · xₖ: accepted numerical state; x₀≠0

The negative sign still alternates, but the magnitude above one now causes growth. $$a=1-4(0.6)=-1.4,\qquad |a|=1.4>1.$$

Solution 5 / 5

Answer · C

h: positive time step · a: Euler multiplier · xₖ: accepted numerical state; x₀≠0

A confuses a negative multiplier with growth.
B assumes physical stability guarantees numerical stability.
C checks the three multiplier magnitudes: zero, point six and one point four.

Step-size restriction; asymptotic decay — [4] FNC §11.3 · Course example

Your turn

Numerical Stability

Every mode must satisfy the bound

$$\dot x=Ax,\quad A=\operatorname{diag}(-2,-8)$$

A: constant state matrix
h: one step used for both modes
x₀: arbitrary initial state

For Euler, which interval gives decay from EVERY start?

0<h≤1/4
0<h<1/4
0<h<1

✓

Correct B. Both strict bounds must hold.

Solution 1 / 5

1 · Require both modes to decay

λ₁,λ₂: modal rates −2 and −8 · h: shared positive step size · aᵢ: Euler multiplier 1+hλᵢ for mode i

An arbitrary initial state can contain either mode, so both multiplier magnitudes must be below one. $$|1+h\lambda_i|<1\quad\Longleftrightarrow\quad-1<1+h\lambda_i<1.$$

Solution 2 / 5

2 · Bound the slow mode

λ₁,λ₂: modal rates −2 and −8 · h: shared positive step size · aᵢ: Euler multiplier 1+hλᵢ for mode i

Subtract one, divide by −2, and reverse both inequalities. $$-1<1-2h<1\ \Longrightarrow\ -2<-2h<0\ \Longrightarrow\ 0<h<1.$$

Solution 3 / 5

3 · Bound the fast mode

λ₁,λ₂: modal rates −2 and −8 · h: shared positive step size · aᵢ: Euler multiplier 1+hλᵢ for mode i

The same calculation with −8 gives the more restrictive bound. $$-1<1-8h<1\ \Longrightarrow\ -2<-8h<0\ \Longrightarrow\ 0<h<\tfrac14.$$

Solution 4 / 5

4 · Intersect and check the endpoint

λ₁,λ₂: modal rates −2 and −8 · h: shared positive step size · aᵢ: Euler multiplier 1+hλᵢ for mode i

The endpoint is excluded: its fast mode alternates with unchanged magnitude. $$(0,1)\cap(0,\tfrac14)=(0,\tfrac14),\qquad h=\tfrac14\Rightarrow a_2=-1.$$

Solution 5 / 5

Answer · B

λ₁,λ₂: modal rates −2 and −8 · h: shared positive step size · aᵢ: Euler multiplier 1+hλᵢ for mode i

A includes a multiplier minus one, so decay fails.
B intersects both strict modal bounds.
C ignores the restrictive fast mode.

Stability of a linear system — [4] FNC §11.3 · Course example

Diagnostic

Numerical Stability

Euler adds energy to an oscillator

$$\begin{gathered}\dot x=Ax,\quad A=\begin{bmatrix}0&-1\\1&0\end{bmatrix}\\|1+ih|=\sqrt{1+h^2}>1\quad(h>0)\end{gathered}$$

±i: eigenvalues of A
h: positive Euler step

A bounded exact orbit can produce an outward numerical spiral.

Absolute stability; harmonic oscillator — [4] FNC §11.3

Recap

Numerical Stability

Model validation and solver verification

$$\text{model: }e^{\lambda h}\qquad\text{Euler: }1+h\lambda$$

λ: continuous rate
h: step size
two factors: exact and numerical evolution over the same interval

Invariant: constant along an exact trajectory. Solver verification checks equations; model validation checks measurements.

Verification; validation — [5] NASA, V&V glossary

Why

Runge–Kutta Methods

Sample the field inside the step

$$\Delta x=\int_{t_k}^{t_k+h}f(x(s),s)\,ds$$

Δx: exact increment from the current exact state
f: state-rate model along the exact trajectory
h: one integration interval

A better slope estimate can improve the same time step.

Runge–Kutta methods — [4] FNC §6.4 · Course diagram

Rule

Runge–Kutta Methods

Construct a trial midpoint

$$\begin{aligned}k_1&=f(x_k,t_k),\quad x_{\rm trial}=x_k+\tfrac h2 k_1\\k_2&=f(x_{\rm trial},t_k+\tfrac h2)\end{aligned}$$

k₁,k₂: initial and trial-midpoint rates
xtrial: temporary state at tₖ+h/2
h k₁/2: estimated change over half a step

Explicit midpoint: second-order Runge–Kutta (RK2). FNC also calls it improved Euler (IE2).

Explicit midpoint method (RK2) — [4] FNC §6.4 · Course diagram

Worked example

Runge–Kutta Methods

Advance with the midpoint rate

$$x_{k+1}=x_k+h k_2$$

xₖ: original current state
k₂: rate evaluated at the temporary midpoint state
h: full step length

ẋ=−2x, x₀=1, h=.2: k₁=−2 → trial=.8 → k₂=−1.6 → x₁=.68.

Explicit midpoint method (RK2) — [4] FNC §6.4 · Course diagram

Your turn

Runge–Kutta Methods

Your turn: halve the midpoint step

$$\dot x=-2x,\quad x(0)=x_0=1,\quad h=0.1$$

xtrial: x₀+(h/2)k₁; temporary midpoint state
k₁,k₂: f(x₀), f(xtrial): rates
h k₂: estimated change over the whole step

Explicit midpoint: which x₁ estimates x(0+h)?

.9: xtrial
.72: xtrial + h k₂
.82: x₀ + h k₂

✓

Correct C. Return to the original state.

Solution 1 / 5

1 · Evaluate the initial rate

k₁,k₂: rates · xtrial: temporary state at t=0.05 · h k₂: estimated full-step change · x₁: next state at t=0.1

Use the original state x₀=1. $$k_1=f(x_0)=-2(1)=-2.$$

Solution 2 / 5

2 · Build the trial midpoint

k₁,k₂: rates · xtrial: temporary state at t=0.05 · h k₂: estimated full-step change · x₁: next state at t=0.1

Advance only half a step using the initial rate. $$x_{\rm trial}=x_0+\frac h2k_1=1+0.05(-2)=0.9.$$

Solution 3 / 5

3 · Evaluate the midpoint rate

k₁,k₂: rates · xtrial: temporary state at t=0.05 · h k₂: estimated full-step change · x₁: next state at t=0.1

The second model call uses the trial state. $$k_2=f(x_{\rm trial})=-2(0.9)=-1.8.$$

Solution 4 / 5

4 · Compute the state at t = 0.1

k₁,k₂: rates · xtrial: temporary state at t=0.05 · h k₂: estimated full-step change · x₁: next state at t=0.1

The full-step change is $hk_2=0.1(-1.8)=-0.18$. Add it to the original state. $$x_1=x_0+hk_2=1-0.18=0.82.$$

Solution 5 / 5

Answer · C

k₁,k₂: rates · xtrial: temporary state at t=0.05 · h k₂: estimated full-step change · x₁: next state at t=0.1

A stops at the temporary trial state.
B adds the final increment to the wrong base state.
C uses the midpoint rate but updates from the original state.

Explicit midpoint method (RK2) — [4] FNC §6.4 · Course diagram

Rule

Runge–Kutta Methods

Classical RK4 evaluates four rates

$$\begin{aligned}k_1&=f(x,t)\\k_2&=f(x+hk_1/2,t+h/2)\\k_3&=f(x+hk_2/2,t+h/2)\\k_4&=f(x+hk_3,t+h)\end{aligned}$$

x,t: current state and time
h: one full step
k₁…k₄: four rates with state/time units

RK4: classical fourth-order Runge–Kutta. Each trial starts from the same x.

Classical fourth-order Runge–Kutta (RK4) — [4] FNC §6.4 · Course diagram

Worked example

Runge–Kutta Methods

Combine the four rates

$$x_{k+1}=x_k+\frac h6(k_1+2k_2+2k_3+k_4)$$

k₁…k₄: RK4 rates
1,2,2,1: weights before division by six
h: full step length

ẋ=−2x, x₀=1, h=.2: rates (−2,−1.6,−1.68,−1.328) give x₁=.6704.

Classical RK4; weighted stage sum — [4] FNC §6.4 · Course diagram

Your turn

Runge–Kutta Methods

Your turn: compare equal model-call budgets

$$\dot x=-2x,\quad x_0=1,\quad h=0.1$$

h: one RK4 step interval
cost: number of calls to f

How many Euler substeps match ONE RK4 step in model calls?

Four, each h/4, over the same interval
One, with step h
Four, each h, over the same interval

✓

Correct A. Match cost AND interval.

Solution 1 / 4

1 · Match the call budget

h: physical interval of one RK4 step; h=0.1 · nE,hE: Euler substep count and substep size · model call: one evaluation of f(x)=−2x · x4,E: state after four Euler substeps

Classical RK4 needs four model calls per step; Euler needs one. $$n_E\times1=4\quad\Longrightarrow\quad n_E=4.$$

Solution 2 / 4

2 · Match the physical interval

h: physical interval of one RK4 step; h=0.1 · nE,hE: Euler substep count and substep size · model call: one evaluation of f(x)=−2x · x4,E: state after four Euler substeps

The four Euler substeps must add up to the same total time h. $$4h_E=h\quad\Longrightarrow\quad h_E=h/4=0.025.$$

Solution 3 / 4

3 · Verify the resulting Euler run

h: physical interval of one RK4 step; h=0.1 · nE,hE: Euler substep count and substep size · model call: one evaluation of f(x)=−2x · x4,E: state after four Euler substeps

For this scalar model, each substep multiplies the state by 0.95. $$x_{4,E}=(1-2h_E)^4x_0=(1-0.05)^4(1)=0.95^4=0.81450625.$$

Solution 4 / 4

Answer · A

h: physical interval of one RK4 step; h=0.1 · nE,hE: Euler substep count and substep size · model call: one evaluation of f(x)=−2x · x4,E: state after four Euler substeps

A matches both cost and interval.
B matches the interval but uses only one call.
C matches cost but advances four times as far.

Function evaluations; order of accuracy — [4] FNC §6.4

Comparison

Runge–Kutta Methods

Accuracy improves at a computational price

$$\varepsilon_{\rm Euler}=O(h),\quad\varepsilon_{\rm midpoint}=O(h^2),\quad\varepsilon_{\rm RK4}=O(h^4)$$

ε: norm of state error at a fixed final time
h: step tending to zero
assumptions: bounded derivatives; field locally Lipschitz near solution

Absolute-stability region: hλ values with bounded repeated updates on ẋ=λx. A large RK4 step can fail.

Order of accuracy; absolute stability — [4] FNC §6.4 · [4] FNC §11.3

Model

Pendulum Simulation

The restoring acceleration has a sign

$$\dot q=v,\qquad\dot v=-\frac{g_0}{\ell}\sin q$$

q: angle from downward vertical
v: angular velocity
g₀,ℓ: positive gravity acceleration and length

Unforced, undamped motion is the model for this comparison.

Simple pendulum; first-order system — [1] Tedrake, ch.2

Definition

Pendulum Simulation

Mechanical energy is a diagnostic

$$E(q,v)=\tfrac12m\ell^2v^2+mg_0\ell(1-\cos q)$$

q,v: angle and angular velocity
m,ℓ,g₀: mass, length and gravity
two terms: kinetic and potential energy

The zero of potential energy is a convention; energy differences have meaning.

Mechanical energy — [1] Tedrake, ch.2

Derivation

Pendulum Simulation

The exact dynamics conserve energy

$$\begin{aligned}E&=\tfrac12m\ell^2v^2+mg_0\ell(1-\cos q)\\\dot E&=m\ell^2v\dot v+mg_0\ell\sin q\,\dot q=0\end{aligned}$$

q̇=v: kinematic identity
v̇: −(g₀/ℓ)sin q
Ė: rate of mechanical energy

The two terms cancel only for the stated unforced, undamped model.

Conservation of mechanical energy — [1] Tedrake, ch.2

Your turn

Integrated Pendulum Problem

Choose the state and the check first

$$\ddot q=-\sin q,\quad q_0=1,\quad v_0=0$$

q₀,v₀: angle and angular velocity
model: unforced, undamped; m=ℓ=g₀=1

Which state and field give a first-order ODE for this model?

x=q; f=−sin q
x=(q,v)ᵀ; f=(v,−sin q)ᵀ
x=(q,v)ᵀ; f=(−sin q,v)ᵀ

✓

Correct B. Keep state/rate order.

Solution 1 / 4

1 · Introduce velocity

q,v: angle and angular velocity · x: complete column state · f(x): state derivative for normalized undamped dynamics

A second-order equation needs both position and velocity to determine the future. $$v=\dot q,\qquad\ddot q=-\sin q\quad\Longrightarrow\quad\dot v=-\sin q.$$

Solution 2 / 4

2 · Stack the state and its rates

q,v: angle and angular velocity · x: complete column state · f(x): state derivative for normalized undamped dynamics

Use the same coordinate order in x and f. $$\begin{gathered}x=\begin{bmatrix}q\\v\end{bmatrix},\\\dot x=f(x)=\begin{bmatrix}\dot q\\\dot v\end{bmatrix}=\begin{bmatrix}v\\-\sin q\end{bmatrix}.\end{gathered}$$

Solution 3 / 4

3 · Apply the given initial state

q,v: angle and angular velocity · x: complete column state · f(x): state derivative for normalized undamped dynamics

The pendulum starts at angle one with zero angular velocity. $$x_0=(1,0)^T,\qquad f(x_0)=(0,-\sin1)^T\approx(0,-0.84147098)^T.$$

Solution 4 / 4

Answer · B

q,v: angle and angular velocity · x: complete column state · f(x): state derivative for normalized undamped dynamics

A omits velocity and mistakes acceleration for angle rate.
B carries both state coordinates and returns their rates in order.
C swaps the rate components.

Pendulum IVP; numerical integration — [1] Tedrake, ch.2 · [4] FNC §6.4 · Course example

Your turn

Integrated Pendulum Problem

Compare states at t = 0.1

$$x_0=(1,0)^T,\ h=.1,\ f(q,v)=(v,-\sin q)^T$$

x=(q,v)ᵀ: angle and angular velocity
h: full time step

Why can midpoint give q₁≠q₀ even though v₀=0?

Nonzero velocity at the trial midpoint.
Twice the physical time.
A different initial velocity.

✓

Correct A. Trial velocity need not be zero.

method	q₁	v₁
Euler	1.000000	−.084147
midpoint	.995793	−.084147
RK4	.995795	−.084071

Solution 1 / 12