Ordinary Differential Equation (ODE)

I’ll help you solve an Ordinary Differential Equation (ODE) using $Python$.

Let’s solve a spring-mass-damper system, which is a classic second-order ODE problem.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def spring_mass_damper(state, t, m, k, c):
"""
Spring-mass-damper system
Second order ODE: m*d²x/dt² + c*dx/dt + k*x = 0
Converted to system of first order ODEs:
dx/dt = v
dv/dt = -(k/m)x - (c/m)v
"""
x, v = state
dxdt = v
dvdt = -(k/m) * x - (c/m) * v
return [dxdt, dvdt]

# System parameters
m = 1.0 # mass (kg)
k = 10.0 # spring constant (N/m)
c = 0.5 # damping coefficient (Ns/m)

# Initial conditions
x0 = 1.0 # initial displacement (m)
v0 = 0.0 # initial velocity (m/s)
state0 = [x0, v0]

# Time points
t = np.linspace(0, 10, 1000)

# Solve ODE
solution = odeint(spring_mass_damper, state0, t, args=(m, k, c))
x = solution[:, 0] # displacement
v = solution[:, 1] # velocity

# Calculate energy
potential_energy = 0.5 * k * x**2
kinetic_energy = 0.5 * m * v**2
total_energy = potential_energy + kinetic_energy

# Create plots
plt.figure(figsize=(15, 10))

# Position vs Time
plt.subplot(2, 2, 1)
plt.plot(t, x, 'b-', label='Position')
plt.grid(True)
plt.xlabel('Time (s)')
plt.ylabel('Position (m)')
plt.title('Position vs Time')
plt.legend()

# Velocity vs Time
plt.subplot(2, 2, 2)
plt.plot(t, v, 'r-', label='Velocity')
plt.grid(True)
plt.xlabel('Time (s)')
plt.ylabel('Velocity (m/s)')
plt.title('Velocity vs Time')
plt.legend()

# Phase Portrait
plt.subplot(2, 2, 3)
plt.plot(x, v, 'g-', label='Phase Portrait')
plt.grid(True)
plt.xlabel('Position (m)')
plt.ylabel('Velocity (m/s)')
plt.title('Phase Portrait')
plt.legend()

# Energy vs Time
plt.subplot(2, 2, 4)
plt.plot(t, potential_energy, 'b-', label='Potential Energy')
plt.plot(t, kinetic_energy, 'r-', label='Kinetic Energy')
plt.plot(t, total_energy, 'g-', label='Total Energy')
plt.grid(True)
plt.xlabel('Time (s)')
plt.ylabel('Energy (J)')
plt.title('Energy vs Time')
plt.legend()

plt.tight_layout()
plt.show()

Let me explain the problem and solution in detail:

  1. The Problem: Spring-Mass-Damper System

    • Mathematical form: $m(d²x/dt²) + c(dx/dt) + kx = 0$
    • Where:
      • $x$ is displacement
      • $m$ is mass
      • $k$ is spring constant
      • $c$ is damping coefficient
  2. Solution Method:

    • Convert 2nd order ODE to system of two 1st order ODEs
    • Use scipy’s odeint for numerical integration
    • System parameters:
      • Mass = $1.0 kg$
      • Spring constant = $10.0 N/m$
      • Damping coefficient = $0.5 Ns/m$
  3. Visualization (four plots):

  • Position vs Time:

    • Shows oscillatory motion with damping
    • Amplitude decreases over time due to damping
  • Velocity vs Time:

    • Shows velocity oscillations
    • Amplitude decreases faster than position
  • Phase Portrait:

    • Shows system trajectory in position-velocity space
    • Spiral pattern indicates damped oscillation
  • Energy vs Time:

    • Shows potential and kinetic energy exchange
    • Total energy decreases due to damping
  1. Physical Interpretation:

    • System starts at $x = 1m$ with zero velocity
    • Oscillates back and forth while losing energy
    • Motion gradually decreases due to damping
    • System eventually comes to rest at equilibrium ($x = 0$)
  2. Energy Analysis:

    • Energy constantly transfers between potential and kinetic
    • Total energy decreases exponentially due to damping
    • Potential energy is maximum at maximum displacement
    • Kinetic energy is maximum at zero displacement