Exact solution for path finding example.

b.py

import numpy as np
import pylab as plt
from skfmm import travel_time, distance
from scipy.interpolate import RectBivariateSpline
from scipy.optimize import fminbound
from math import sin

# Brachistochrone curve
# see: http://mathworld.wolfram.com/BrachistochroneProblem.html
b=1/np.pi
theta = np.linspace(0,np.pi,300)
x = b*(theta-np.sin(theta))
y = b*(1-np.cos(theta))


N     = 111 # value can be 111 or 1112
assert (np.linspace(0,1.1,N)==1.0).any()
X, Y  = np.meshgrid(np.linspace(0, 1.1, N), np.linspace(0, 1.1, N))
grid_spacing = 1.1/float(N)
phi = np.ones_like(X)
phi[X==1] = 0
phi[X>1] = -1
dist = distance(phi)

vel = np.sqrt(Y)
time = travel_time(phi,vel,dx=grid_spacing)
plt.contourf(X,Y,time)
plt.plot(x,y) # plot exact solution as solid line

# starting point for numerical solution
# we cannot start a the origin because the velocity is zero.
y0 = 0.05
theta0 = np.arccos(1-y0/b)
x0 = b*(theta0-np.sin(theta0))
print y0, x0
plt.plot([x0], [y0], "*") # plot starting point

def optimal_path_2d(travel_time, starting_point, dx, coords,
                 N=100000):
    """
    Find the optimal path from starting_point to the zero contour
    of travel_time. dx is the grid spacing
    Solve the equation x_t = - grad t / | grad t |
    """
    grad_t_y, grad_t_x = np.gradient(travel_time, dx)
    if isinstance(travel_time, np.ma.MaskedArray):
        grad_t_y[grad_t_y.mask]=0.0
        grad_t_y = grad_t_y.data
        grad_t_x[grad_t_x.mask]=0.0
        grad_t_x = grad_t_x.data

    gradx_interp = RectBivariateSpline(coords, coords,
                                       grad_t_x)
    grady_interp = RectBivariateSpline(coords, coords,
                                       grad_t_y)

    def get_velocity(position):
        """ return normalized velocity at pos """
        x, y = position
        vel = np.array([gradx_interp(y, x)[0][0],
                        grady_interp(y, x)[0][0]])
        return vel / np.linalg.norm(vel)

    def euler_point_update(pos, ds):
        return pos - get_velocity(pos) * ds

    def runge_kutta(pos, ds):
        """ Fourth order Runge Kutta point update """
        k1 = ds * get_velocity(pos)
        k2 = ds * get_velocity(pos - k1/2.0)
        k3 = ds * get_velocity(pos - k2/2.0)
        k4 = ds * get_velocity(pos - k3)
        return pos - (k1 + 2*k2 + 2*k3 + k4)/6.0

    x = runge_kutta(starting_point, dx)
    xl, yl = [x[1]], [x[0]]
    for i in range(N):
        xl.append(x[1])
        yl.append(x[0])
        x = runge_kutta(x, dx)
        distance = ((x[1] - xl[-1])**2 + (x[0] - yl[-1])**2)**0.5
        if distance < dx * 0.9999:
            print "exiting"
            break
        # We should come up with a better stopping criteria that puts
        # a point exactly on the zero contour.

    return yl, xl

coords = np.linspace(0, 1.1, N)
dx = coords[1]-coords[0]
xp, yp = optimal_path_2d(time, (x0,y0), dx, coords)
plt.plot(xp, yp, "o") # plot numerical solution as points.
plt.gca().set_aspect(1)
plt.axvline(1)
plt.colorbar()


# this is one way to find the error but it is problematic near the end
# point in this case.
thetap = np.arccos(1-np.array(yp)/b)
x_exact = b*(thetap-np.sin(thetap))
error = xp-x_exact
print "error metric 1"
print error

def find_theta(x):
    func = lambda theta : abs(x-b*(theta-sin(theta)))
    return fminbound(func, 0, np.pi, xtol=1e-12)

thetap2 = np.array(map(find_theta, xp))
# this number should be smaller than the errors calculated
print abs(xp - b*(thetap2-np.sin(thetap2))).max()
y_exact = b*(1-np.cos(thetap2))
print "error metric 2"
print yp-y_exact

print "==== analytic solution to travel time as a function of theta"
# first find the arc-length of the cycloid as a function of theta
# http://www-math.mit.edu/~djk/18_01/chapter18/section02.html
arc_length = 4/np.pi*(1-np.cos(np.pi/2))
dx=x[1:]-x[:-1]
dy=y[1:]-y[:-1]


print "arc-length exact, numerical", arc_length, sum(np.sqrt(dx**2+dy**2))

# next find travel time along this arc-length
# integrate(2*(sin(a/2))/(sqrt((1-cos(a))/pi))/pi,a)
# this is travel time as a function of theta
tt= lambda a : 2*a*np.sin(a/2.0)/np.sqrt(np.pi-np.pi*np.cos(a))
tt0 = tt(np.pi) # total travel time along Brachistochrone curve

ybar=(y[1:]+y[:-1])/2.0
ntt=sum(np.sqrt(dx**2+dy**2)/np.sqrt(ybar))
print "travel time exact, numerical", tt0, ntt

tstar = tt0-tt(np.pi/2.0)
xx = b*(np.pi/2-np.sin(np.pi/2))
yy = b*(1-np.cos(np.pi/2))
print "travel time at {}, {} : {}".format(xx,yy,tstar)
plt.plot((xx,),(yy,), "o")

c=0.3
b=(1-c)/np.pi
x = b*(theta-np.sin(theta))+c
y = b*(1-np.cos(theta))
plt.plot(x,y,linewidth=3)


y0 = 0.05
theta0 = np.arccos(1-y0/b)
x0 = b*(theta0-np.sin(theta0))+c
print y0, x0
plt.plot([x0], [y0], "*") # plot starting point
dx = coords[1]-coords[0]
xp, yp = optimal_path_2d(time, (x0,y0), dx, coords)
plt.plot(xp, yp, "o") # plot numerical solution as points.


plt.show()



# how to use this to check the accuracy of travel_time???

# ttt = tt(theta)
# plt.plot(theta, ttt)
# plt.show()

import numpy as np
import pylab as plt

from math import pi, acos, sin
def residual(a):
    thetaa = acos(1-yb*pi/(1-a))
    return abs(abs((1-a)/pi*(thetaa-sin(thetaa))+a)-abs(xb))

yb = 0.4
xb = 0.4
from scipy.optimize import fminbound
ret = fminbound(residual,0,0.3,full_output=True, xtol=1e-12)
a = ret[0]
print a
#theta = acos(1-yb*pi/(1-a))
#print theta

plt.plot((xb,),(yb,),"o")
theta = np.linspace(0,pi,2000)
c=a
b=(1-c)/np.pi
x = b*(theta-np.sin(theta))+c
y = b*(1-np.cos(theta))
plt.gca().set_aspect(1)
plt.plot(x,y,"o-")
plt.show()

# yhat = (1-xhat)*2/pi