Uwzględnienie niezerowej siły wymuszającej w kodzie programu jest nieskomplikowane – wystarczy uzupełnić linię 11 kodu, uzupełniając funkcję o dodatkowy parametr – aktualną wartość siły wymuszającej.
Dodatkowo program generuje cztery wykresy zależności gdzie – współczynnik tłumienia, korzystając online z biblioteki mathplotlib.pyplot.
import math import matplotlib.pyplot as plt
g = 9.81 # Przyspieszenie ziemskie L = 9.81 # Długość wahadła
T_0 = 2 * math.pi * math.sqrt(L / g)
def f(theta, omega, c): # Równania ruchu dla wahadła z tłumieniem omega_dot = (-c * omega - (g / L) * math.sin(theta)) return omega_dot
def runge_kutta(theta, omega, h, c): k1_theta = h * omega k1_omega = h * f(theta, omega, c)
k2_theta = h * (omega + k1_omega / 2) k2_omega = h * f(theta + k1_theta / 2, omega + k1_omega / 2, c)
k3_theta = h * (omega + k2_omega / 2) k3_omega = h * f(theta + k2_theta / 2, omega + k2_omega / 2, c)
k4_theta = h * (omega + k3_omega) k4_omega = h * f(theta + k3_theta, omega + k3_omega, c)
theta += (k1_theta + 2 * k2_theta + 2 * k3_theta + k4_theta) / 6 omega += (k1_omega + 2 * k2_omega + 2 * k3_omega + k4_omega) / 6
return theta, omega
def main(): # DANE: theta_0 = 170 # Początkowy kąt wychylenia [w stopniach]
h = 0.00001 # Krok czasu iterations = 4000000 # Liczba iteracji