Arc length

2003 days ago by JAVIER

import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from scipy import integrate 
       
N = 10 a = 0 b = np.pi var('x') #---------------------- fun = cos(pi*x) #fun = x*log(x+1)*cos(x^2/pi) #---------------------- length_fun= sqrt(1 + fun.diff(x)^2) #---------------------- fun_f = fast_callable(fun,vars=[x]) length_fun_f = fast_callable(length_fun,vars=[x]) #---------------------- exact_length = integrate.quad(length_fun_f,a,b) print 'Exact length : ', exact_length[0] #---------------------- x = np.linspace(a,b,100) y = fun_f(x) fig = plt.figure() axes=fig.add_subplot(111) axes.plot(x,y,'-',lw = 2) #----------------------partition xp = np.linspace(a,b,N,endpoint=True) yp = fun_f(xp) axes.plot(xp,yp,'ko') L = [] #---------------------- for h,k in zip(xp,yp): L.append((h,k)) #---------------------- for k in range(len(L)-1): plt.plot([L[k][0], L[k+1][0]], [L[k][1], L[k+1][1]], 'k-', lw=2) #---------------------- fig.set_size_inches(4,4) fig.savefig('test.png') plt.close() #---------------------- def distance(p1,p2): return sqrt((p1[0]-p2[0])^2 + (p1[1]-p2[1])^2) #---------------------- approx_length = 0 for k in range(len(L)-1): approx_length += distance(L[k],L[k+1]) print 'Approx. length: ', approx_length.n(20) 
       
Exact length  :  7.09272996122
Approx. length:  6.7946
Exact length  :  7.09272996122
Approx. length:  6.7946