generating a nonuniform mesh

mesh.py

#!/usr/bin/env python

from scipy import optimize
from scipy import integrate
import numpy as np
import sys

l, L = 0.02, 1
N = 128

mu, sigma = 0.5, 0.1

# density distrubtion function
def rho(x):
    return 1.0 + np.tanh((x-(mu-sigma))/l) - np.tanh((x-(mu+sigma))/l)

# indefinite integral of the density distrubtion function
def rhoint(x):
    return x + l*np.log(np.cosh((x - mu + sigma)/l)) - l*np.log(np.cosh((-x+mu+sigma)/l))


# first print density distribution function
for x in np.linspace(0, 1, N):
    print x, rho(x)
print

# and find normalisation constant
theta = integrate.quad(rho, 0, L)[0]
print >> sys.stderr, "theta =", theta




def f(x, x0):
    return rhoint(x[0]) - rhoint(x0) - theta/N

xs = np.zeros(N)

x0 = 0
for i in xrange(N):
    x0 = optimize.fsolve(f, x0, args=(x0,))[0]
    xs[i] = x0

    print x0, np.sin(np.pi * x0)