Matlab ode45 tutorial and exercises

This tutorial covers the implementation of a mass-spring-damper system and various other physical systems using ode45. Use the verification tools to match your graphs for all state variables.

Part 1: Guided Tutorial (Mass-Spring-Damper)

Mathematical Background

The 2nd-order ODE for a mass-spring-damper system is $m\ddot{x} + c\dot{x} + kx = 0$. Using state-space variables $y_1 = x$ and $y_2 = \dot{x}$:

$$\dot{y}_1 = y_2$$ $$\dot{y}_2 = \frac{1}{m}(-cy_2 - ky_1)$$

MATLAB Implementation

First, define the dynamics in msd_dynamics.m:

function dydt = msd_dynamics(t, y, m, c, k)

    % y(1) = position, y(2) = velocity

    dydt = [y(2); (1/m)*(-c*y(2) - k*y(1))];

end

Next, create the Runner Script run_msd.m to solve and plot all states:

% run_msd.m - Simulation Entry Point

m = 1.0; c = 0.5; k = 10.0; 

tspan = [0 15]; 

y0 = [1.0; 0.0]; % [Initial Position; Initial Velocity]

[t, y] = ode45(@(t,y) msd_dynamics(t, y, m, c, k), tspan, y0);

% Plotting all states (Time History)

plot(t, y(:,1), 'b-', t, y(:,2), 'r--');

legend('Position (m)', 'Velocity (m/s)');

xlabel('Time (s)'); ylabel('Amplitude'); grid on;

To execute the simulation, type this into the MATLAB Command Window:

run_msd

Mass (m)

Damping (c)

Spring (k)

Time (s)	Pos (m)	Vel (m/s)

Exercise 1: Van der Pol Oscillator

Solve for $\ddot{x} - \mu(1 - x^2)\dot{x} + x = 0$. Provide time histories for both $x$ and $\dot{x}$.

Nonlinearity ($\mu$)

Init $x_0$

Time	$y_1$ (Pos)	$y_2$ (Vel)

Exercise 2: Non-Linear Pendulum

Model the pendulum $\ddot{\theta} + \frac{g}{L} \sin(\theta) = 0$. Plot angle $\theta$ and angular velocity $\omega$.

Length (L)

Angle ($\theta_0$)

Time	$\theta$ (rad)	$\omega$ (rad/s)

Exercise 3: Duffing Oscillator

Model $\ddot{x} + \delta\dot{x} + \beta x + \alpha x^3 = 0$. Display both states over 30 seconds.

Damping ($\delta$)

Cubic ($\alpha$)

Time	Pos	Vel

Exercise 4: Lorenz System (3D Chaos)

This requires a 3-element state vector: $y = [x, y, z]^T$. Plot all three time histories.

Contextual Overview: The Lorenz system models atmospheric convection. It is famously sensitive to initial conditions.

$\sigma$: The Prandtl number.
$\rho$: The Rayleigh number (determines the "chaos" level).
State Variables: $x$ (convection rate), $y$ (temp difference), $z$ (distortion of temp profile).

$$\dot{x} = \sigma(y - x)$$ $$\dot{y} = x(\rho - z) - y$$ $$\dot{z} = xy - \beta z$$

Sigma ($\sigma$)

Rho ($\rho$)

Time	X	Y	Z

Exercise 5: Restricted 3-Body Problem (CR3BP)

Mathematical & Physical Context

The Circular Restricted Three-Body Problem (CR3BP) models the motion of a negligible mass (satellite) under the influence of two large primary masses (e.g., Earth and Moon) that orbit each other. We use a rotating reference frame where the two primaries appear stationary on the x-axis.

Parameter and State Definitions:

$\mu$ (Mass Ratio): Defined as $m_2 / (m_1 + m_2)$. For the Earth-Moon system, $\mu \approx 0.01215$.
States $y_1, y_2$: Represent the $(x, y)$ position in the rotating plane. The primary mass is at $(-\mu, 0)$ and the secondary is at $(1-\mu, 0)$.
States $y_3, y_4$: Represent the velocity components $(\dot{x}, \dot{y})$ relative to the rotating frame.
Coriolis Acceleration: The terms $2\dot{y}$ and $-2\dot{x}$ are "apparent forces" required because our coordinate system is spinning.
Distances $r_1, r_2$: These are the scalar distances from the satellite to the primaries.
- $r_1 = \sqrt{(x + \mu)^2 + y^2}$ (Distance to the larger primary)
- $r_2 = \sqrt{(x - 1 + \mu)^2 + y^2}$ (Distance to the secondary primary)

$$\ddot{x} = x + 2\dot{y} - \frac{(1-\mu)(x+\mu)}{r_1^3} - \frac{\mu(x-1+\mu)}{r_2^3}$$ $$\ddot{y} = y - 2\dot{x} - \frac{(1-\mu)y}{r_1^3} - \frac{\mu y}{r_2^3}$$

Mass Ratio ($\mu$)

Init $x_0$

Init $y_0$

Init $\dot{x}_0$

Init $\dot{y}_0$

Time	X	Y	$\dot{x}$	$\dot{y}$

Test Scenarios & Expected Results:

1. The $L_4$ Equilibrium (Default): Set $x_0 = 0.48785$, $y_0 = 0.86603$, and all velocities to $0$. Result: The satellite should remain perfectly stationary in the rotating frame (flat lines in time-history).
2. Horseshoe Orbit: Use the Earth-Moon $\mu = 0.01215$. Set $x_0 = 0.99$, $y_0 = 0$, $\dot{x}_0 = 0$, $\dot{y}_0 = -0.01$. Result: The $x$ vs $y$ plot will show a large, "C" shaped path that avoids the secondary mass, characteristic of co-orbital objects.
3. Lyapunov Orbit (Small Oscillation): Set $x_0 = 0.8$, $y_0 = 0$, $\dot{x}_0 = 0$, $\dot{y}_0 = 0.2$. Result: You should see a closed periodic loop around the $L_1$ or $L_2$ region.
4. Chaotic Escape: Increase $\dot{y}_0$ slightly from the Lyapunov case (e.g., $\dot{y}_0 = 0.5$). Result: The satellite will likely "sling-shot" around a primary and be ejected from the system.