import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
from numba import jit, njit
# Konfiguriere matplotlib um schönere Plots zu erhalten
plt.rcParams.update({
"interactive": True,
"figure.dpi": 200,
"figure.figsize": (10, 6),
"font.family": "serif",
"mathtext.fontset": "cm",
"axes.grid": True,
"axes.labelsize": 14,
"grid.linewidth": .5,
"grid.alpha": .5,
"lines.linewidth": 1,
"lines.dashed_pattern": (6, 4),
})
def simulate(f, y0, duration, steps=1000, **params):
h = duration / steps
t_arr = np.linspace(0, h * steps, steps)
y_arr = np.empty((steps, *np.shape(y0)))
y_arr[0] = y0
for i in range(steps-1):
t, y = t_arr[i], y_arr[i]
k1 = f(t, y, **params)
k2 = f(t + h/2, y + h/2 * k1, **params)
k3 = f(t + h/2, y + h/2 * k2, **params)
k4 = f(t + h, y + h * k3, **params)
y_arr[i+1] = y + h/6 * (k1 + 2*k2 + 2*k3 + k4)
return t_arr, y_arr
def spring_f(t, y, alpha=1, beta=.05, gamma=1, delta=.1, omega=1):
"""Creates a derivative function describing a special spring for the given params."""
return np.array([
y[1],
gamma * np.cos(omega * t) - delta * y[1] - alpha * y[0] - beta * y[0]**3
])
def plot_spring(t, y):
x, v = y[:,0], y[:,1]
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, dpi=200)
ax1.plot(t, x, c="b")
ax1.set_ylabel("$x$")
ax2.plot(t, v, c="r")
ax2.set_ylabel("$v$")
ax2.set_xlabel("$t$")
t, y = simulate(spring_f, y0=np.array([0, 0]), duration=100)
plot_spring(t, y)
def simulate_spring_periods(omega, periods, steps_per_period, y0=np.zeros(2), **params):
duration = periods * 2 * np.pi / omega
steps = periods * steps_per_period
return simulate(spring_f, y0, duration, steps, omega=omega, **params)
def determine_amps_and_phases(omegas, steps_per_period=100, periods=2, start_periods=10, **params):
# Einschwingphase simulieren
start_t, start_y = simulate_spring_periods(omegas[0], start_periods, steps_per_period, **params)
# Ampltituden und Phasen für alle omegas messen
y0 = start_y[-1]
amps = np.empty(len(omegas))
phases = np.empty(len(omegas))
for i, omega in enumerate(omegas):
t, y = simulate_spring_periods(omega, periods, steps_per_period, y0=y0, **params)
# Set x0 and v0 for next omega simulation
y0 = y[-1]
x = y[:,0]
amps[i] = np.max(x)
intercept_index = np.argmax(x)
intercept_t = t[intercept_index]
intercept_phase = (intercept_t * omega) % (2 * np.pi)
phases[i] = intercept_phase
return amps, phases
def plot_amps_and_phases(omegas, amps, phases):
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, dpi=200)
ax1.scatter(omegas, amps, s=1, c="r")
ax1.set_ylabel("$A$")
ax2.scatter(omegas, phases, s=1, c="b")
ax2.set_ylabel("$\\phi$")
ax2.set_xlabel("$\\omega$")
omegas = np.linspace(1/2, 2, 100)
amps, phases = determine_amps_and_phases(omegas)
plot_amps_and_phases(omegas, amps, phases)
def analyze_betas():
betas = [0, 0.01, 0.04, 0.1]
#betas = [.04]
colors = ["r", "b", "g", "y"]
omegas = np.linspace(1/2, 2, 200)
rev_omegas = omegas[::-1]
fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, dpi=200)
fig.set_size_inches(12, 8)
ax1.set_ylabel("$A$")
ax2.set_ylabel("$\\phi$")
ax2.set_xlabel("$\\omega$")
for beta, color in zip(betas, colors):
forward_amps, forward_phases = determine_amps_and_phases(omegas, beta=beta)
ax1.scatter(omegas, forward_amps, c=color, s=8, marker=">", label=f"$\\beta = {beta}$ (vorwärts)")
ax2.scatter(omegas, forward_phases, c=color, s=8, marker=">")
backward_amps, backward_phases = determine_amps_and_phases(rev_omegas, beta=beta)
ax1.scatter(rev_omegas, backward_amps, c=color, alpha=.5, s=8, marker="<", label=f"$\\beta = {beta}$ (rückwärts)")
ax2.scatter(rev_omegas, backward_phases, c=color, alpha=.5, s=8, marker="<")
ax1.legend()
analyze_betas()
def analyze_beta_dependency():
betas = np.linspace(0, 0.1, 50)
omegas = np.linspace(1/2, 2, 200)
rev_omegas = omegas[::-1]
forward_res = np.empty_like(betas)
backward_res = np.empty_like(betas)
for i, beta in enumerate(betas):
forward_amps, forward_phases = determine_amps_and_phases(omegas, steps_per_period=20, beta=beta)
backward_amps, backward_phases = determine_amps_and_phases(rev_omegas, steps_per_period=20, beta=beta)
forward_res[i] = omegas[np.argmax(forward_amps)]
backward_res[i] = rev_omegas[np.argmax(backward_amps)]
fig, ax = plt.subplots()
fig.set_size_inches(10, 6)
ax.scatter(betas, forward_res, s=5, label="vorwärts")
ax.scatter(betas, backward_res, s=5, label="rückwärts")
ax.legend()
analyze_beta_dependency()