Tutorial A12 – Solutions

1 How does logging behave when you rerun the script without restarting the kernel? What happens if you delete the log-file before you rerun?

If you rerun your script, the new logging messages are append to the log file If you delete the log-file in between, the result is the same. The logger remenbers everything written to the log file before and the teh log file is complete restored

2 Play around with the datefmt keyword when configuring the logger to omit displaying year and months.

logging.basicConfig(
    filename="file",
    datefmt='%d-%I:%M:%S %p',
    )

3 Below you find a script that propagates a particle in a harmonic potential numerically (using the euler integrator). You do not need to know about details of the scientific background yet. Just note so much: The script monitors the total energy of the particle at every propagation step. Whenever the change in energy is larger than 0.3, the velocity should be rescaled (multiplied) by a factor of 0.9. This part is not yet included in the script below. If this happens a message should be written to a log file, stating the old and the new speed. Your task is it, to setup a logger and insert a logging statement at the right position in the code.

import logging

fname = "ea8/log.log"
logging.basicConfig(
    filename=fname,
    filemode="w",
    datefmt="%H:%M:%S",
    level=logging.WARNING,
    )


def force(x):
    """Compute the instantaneous force acting on the particle

    Args:
        x: Current position
        x0 (global): Position of the minimum of the potential
        k (global): Force constant
    """

    return -k * (x - x0)


def e_pot(x):
    """Compute the potential energy of the particle

    Args:
        x: Current position
        x0 (global): Position of the minimum of the potential
        k (global): Force constant
    """

    return 0.5 * k * (x - x0)**2


def e_kin(v):
    """Compute the kinetic energy of the particle

    Args:
        v: Particle velocity
        m (global): Particle mass
    """

    return 0.5 * m * v**2


def propagate(x, v):
    """Compute the position of the particle one timestep further

    Args:
        x: Current position
        v: Current velocity
        dt (global): Propagation timestep
        m (global): Particle mass

    Returns:
        New position and velocity
    """

    x_new = x + v*dt + force(x)/(2*m) * dt**2
    v_new = v + force(x)/m * dt

    return x_new, v_new


# Preparing lists to store the computed trajectories
x_t = []  # Position
v_t = []  # Velocity
e_t = []  # Energy

# Particle/Simulation parameters
x0 = 0  # Position of the potential minimum
k = 1.  # Force constant
m = 1.  # Particle mass

dt = 1e-1  # Simulation timestep
T = 1000  # Total number of steps
scale_factor = 0.9
# Rescaling factor to counter energy conservation problems

x = 1  # starting position
v = -2  # starting velocity

# Put them in the output lists
x_t.append(x)
v_t.append(v)
e_t.append(e_pot(x) + e_kin(v))

# Simulate
for i in range(T):
    # Update position and velocity
    x, v = propagate(x, v)
    x_t.append(x)
    v_t.append(v)

    # Update total energy
    e_t.append(e_pot(x) + e_kin(v))

    if (e_t[-1] - e_t[-2]) > 0.3:
        logging.warning(
            'Energy not conserved!'
            f'Rescaling velocity from {v_t[-1]} to {v_t[-1] * scale_factor}'
        )

        v *= scale_factor

# Print the first 10 lines of the log file
with open(fname) as f:
    for line in f.readlines()[:10]:
        print(line, end="")

WARNING:root:Energy not conserved!Rescaling velocity from -7.952943460471916 to -7.157649114424725
WARNING:root:Energy not conserved!Rescaling velocity from 7.818357041905684 to 7.036521337715116
WARNING:root:Energy not conserved!Rescaling velocity from 7.937734070995736 to 7.143960663896163
WARNING:root:Energy not conserved!Rescaling velocity from -7.991841670210238 to -7.192657503189214
WARNING:root:Energy not conserved!Rescaling velocity from 7.831076595508824 to 7.047968935957941
WARNING:root:Energy not conserved!Rescaling velocity from 7.945001174756476 to 7.150501057280829
WARNING:root:Energy not conserved!Rescaling velocity from -8.02546516176625 to -7.222918645589626
WARNING:root:Energy not conserved!Rescaling velocity from 7.8738492791246015 to 7.0864643512121415
WARNING:root:Energy not conserved!Rescaling velocity from 7.9490910738210205 to 7.154181966438919
WARNING:root:Energy not conserved!Rescaling velocity from -8.06551535987205 to -7.258963823884845

4 Plot the result of the simulation as velocity vs. position.

import matplotlib.pyplot as plt

plt.plot(x_t, v_t)
plt.xlabel("$x$")
plt.ylabel("$v$")

Text(0, 0.5, '$v$')