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()