What is the average distance between two randomly chosen points on a Sierpiński triangle?

Question

Since a Sierpiński triangle is made up of three smaller versions of itself, an algorithm made to pick random points on the triangle might go like this:

1. Select one of the three Sierpiński triangles inside the larger Sierpiński triangle
2. Repeat step 2 with the smaller Sierpiński triangle

It will continue getting more specific and will eventually narrow down to a point. If you continue picking two points like this, what will be the average distance? (Euclidian Metric)

Answer 1

I don't know how to do this analytically, but I did a simulation, assuming that the original triangle is equilateral of side $1$. With one million trials, I got a $99\%$ confidence interval of $$(0.42238756083753565, 0.4234628881174311)$$

I started with $$(1,0,0), (0,1,0), (0,0,1)$$ as the vertices of the triangle in barycentric coordinates. Then $30$ times, I chose one of the coordinates uniformly at random, and computed the vertices of the smaller triangle as the average of the chosen vertex and the old vertices. For example, if the first vertex chosen were $(1,0,0)$ the new triangle would have coordinates $$(1,0,0),\left(\frac12,\frac12,0\right),\left(\frac12,0,\frac12\right)$$ After $30$ samples, I took the average of the vertices as the chosen point.

For each trial, I chose two points, and computed the Euclidean distance between them using the formula from Wikipedia. Since I assumed that the side lengths of the triangle are all $1$, this is simply$$d=\sqrt{\frac{x^2+y^2+z^2}2}$$ where $x,y,z$ are the coordinates of the difference between the two chosen points, treated as vectors.

It seemed to me at first that one should be able to do this analytically. If the coordinates of the $n$th triangle are $(X_n, Y_n, Z_n)$ then $X_{n+1}$ is equally likely to be any of $X_n,\ \frac{X_n+Y_n}2,\ \frac{X_n+Z_n}2$. I haven't made any progress, though. This may just be because I'm not thinking clearly, but as I'm typing this it occurs to me that, since the Sierpiński gasket is a fractal, perhaps the distributions functions of X,Y,Z are singular, so that we can't hope for analytic expressions. I'm going to compute the distribution of $X_n$ for small $n$ and see if that gives me any ideas.

For completeness, here's the simulation script:

from sys import argv
from random import choices
from math import sqrt

def main(trials, samples):
    total = 0
    total2 = 0
    for _  in range(trials):
        x = pickPoint(samples)
        y = pickPoint(samples)
        # Triangle side length are all 1
        v =  [(z[0]-z[1])**2 for z in zip(x,y)]
        d = sqrt(sum(v)/2)  # Euclidean distance
        total += d
        total2 += d*d
    mean = total/trials
    variance = total2/trials - mean**2
    sigma = sqrt(variance/trials)
    return mean, sigma

def pickPoint(samples):
    V = [(1,0,0), (0,1,0), (0,0,1)]
    s = choices(range(3), k=samples)
    for k in s:
        v = V[k]
        V = [tuple((z[0]+z[1])/2 for z in zip(v,V[k])) for k in range(3)]  
    return [(V[0][k]+V[1][k]+V[2][k])/3 for k in range(3)]
    
if __name__ == '__main__':
    trials = int(argv[1])
    samples = int(argv[2])
    mu, sigma = main(trials, samples)
    print(f'Mean = {mu}')
    print(f'sigma = {sigma}')
    print(f'{trials} trials')
    print(f'Sample depth = {samples}')
    delta = 2.576*sigma
    print(f'99% confidence interval: ({mu-delta}, {mu+delta})')

When I ran this with parameters 1000000 30 it produced the output

Mean = 0.4229252244774834
sigma = 0.00020872035712257242
1000000 trials
Sample depth = 30
99% confidence interval: (0.42238756083753565, 0.4234628881174311)