generating a nonuniform mesh
#!/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) | |
# 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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment