The Barnsley Fern is a beautiful fractal that can easily be generated in Python.

If we zoom in on one branch, we see that the pattern is repeated:

The python code follows.

from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt


# Main body of program
points = 10000  # The number of points to use.

X = []  # A list of x-coordinates
Y = []  # A list of y-coordinates

# Set the starting point
point = [0.5, 0.0]
X.append(point[0])
Y.append(point[1])


def new_point(p):
    """
    This function takes in a point (x, y) and generates a new point according
    to the given equations.
    """
    r = np.random.uniform(0, 1)
    if r < 0.02:
        p = [0.5, 0.27*p[1]]
    elif 0.02 <= r <= 0.17:
        p = [-0.139*p[0] + 0.263*p[1] + 0.57, 0.246*p[0] + 0.224*p[1] - 0.036]
    elif 0.17 < r <= 0.3:
        p = [0.17*p[0] - 0.215*p[1] + 0.408, 0.222*p[0] + 0.176*p[1] + 0.0893]
    elif 0.3 < r < 1.0:
        p = [0.781*p[0] + 0.034*p[1] + 0.1075, -0.032*p[0] + 0.739*p[1] + 0.27]

    return p

# Generate a large number of points
for i in range(points):
    point = new_point(point)
    X.append(point[0])
    Y.append(point[1])

# Plot the results
plt.scatter(X, Y, c='g', s=.05)
plt.axis('Off')
plt.axes().set_aspect('equal')
plt.title("Barnsley Fern")
plt.savefig("barnsley_fern.png")
plt.show()

This program computes all the prime numbers up to 10,000 using an efficient algorithm.

Instead of checking if a number n is prime by dividing by all previous numbers, or even all previous primes, this program only divides by all previous primes less than or equal to the square root of n.

from __future__ import division, print_function

primes = [2]
print(2)

for n in range(3, 10001):
    i, k = 0, 0
    while primes[i] <= n**(1/2):    # No need to check against ALL smaller primes
        if n % primes[i] == 0:      # Check if divisible by previous primes
            k = 1
        i += 1
    if k == 0:                      # If good, at it to the list and print it
        primes.append(n)
        print(n)

This simple program shows an animation of a two-pendulum version of Newton’s cradle. One pendulum begins at some starting angle and the other pendulum is stationary. When the moving pendulum slams into the stationary one, its momentum is transferred to the other pendulum.

from visual import *
import numpy as np


# Constants
g = 9.80            # (m/s^2)
L = 10              # Length of the pendulums (m)
initialAngle = 1.2  # In radians


# Create the pendulum bob and rod
pend = sphere(pos=(L*np.sin(initialAngle), -L*np.cos(initialAngle), 0), radius=1, color=color.yellow)
rod = cylinder(pos=(0, 0, 0), axis=(pend.pos[0], pend.pos[1], 0), radius=0.1)
pend2 = sphere(pos=(-2, -L, 0), radius=1, color=color.red)
rod2 = cylinder(pos=(-2, 0, 0), axis=(pend2.pos[0]+2, pend2.pos[1], 0), radius=0.1)


def position(right, t):
    """
    Only one of the pendulums is in motion at a given time. This function
    moves the moving pendulum to its new position. We use the equation:
        theta(t) = theta_0*cos(sqrt(g/L)*t)
    """
    theta = initialAngle*np.cos((g/L)**(1/2)*t)

    if right:
        pend.pos = [L * np.sin(theta), -L * np.cos(theta), 0]  # Update position of bob
        rod.axis = [pend.pos[0], pend.pos[1], 0]  # Update rod's position
    else:
        pend2.pos = [L * np.sin(theta) - 2, -L * np.cos(theta), 0]  # Update position of bob
        rod2.axis = [pend2.pos[0] + 2, pend2.pos[1], 0]  # Update rod's position

    # Once the moving pendulum reaches theta = 0, switch to the other one
    if theta <= 0:
        return False  # Return
    else:
        return True

# Increment time
i = 0
right = True  # The right pendulum is the first in motion
while True:
    rate(200)
    right = position(right, i)
    i += 0.01

If you want to include damping, you can change the line

theta = initialAngle*np.cos((g/L)**(1/2)*t)

to

theta = initialAngle*exp(-mu*t)*np.cos((g/L)**(1/2)*t)

with some value for mu.