float lens = length(@P);
@P.y += ( pow(0.5,lens) - 1.0f) / (0.5-1);

import math# Trapezoidal
# S = 1/2(y0+ 2y1 + 2y2 + 2y3+...+ 2yn-1 + yn)def Trapezoidal(down,up,n,func):if up==down:return 0.0h = float(up-down) / float(n)start = func(down)end = func(up)process = 0.0for dt in xrange(0,n+1,1):if dt == 0 or dt == n:continueprocess += 2 * func(down + dt * h)sum =  (start + end + process) * (h/2.0)return sum# Simpson
# S = h/3(y0 + 4y1 + 2y2 + 4y3 + 2y4 + ... + 2yn-1 + yn)
# func is f(x)
def Simpson(down,up,n,func):if up==down:return 0.0h = float(up-down) / float(n)start = func(down)end = func(up)process = 0.0for dt in xrange(0,n+1,1):if dt == 0 or dt == n:continue# select the 1 3 5 7 9... indexif dt%2 == 1:process += 4 * func(down + dt * h)# select the 2 4 6 8 10... indexif dt%2 == 0:process += 2 * func(down + dt * h)sum =  (start + end + process) * (h/3.0)return sum

import numpy as npdef Euler(func, X0, Y0, perspectiveX, dx = 0.1,debug=True):if perspectiveX == X0:return Y0# improve Euler method# Xn = Xn-1 + dx# Zn = Yn-1 + f(Xn-1, Yn-1) * dx# Yn = Yn-1 + [ f(Xn-1, Yn-1) + f(Xn,Zn) / 2 ] * dx
Xn = X0Zn = Y0Yn = Y0for it in np.arange(X0, perspectiveX, dx):if debug:print "-----------loop------------------" , itXn = Xn + dx # 1:X0 + dx , 2:X0 + dx + dx  ...if debug:print "Xn" , XnZn = Yn + func(Yn) * dxYn = Yn + ( (func(Yn) + func(Zn)) / 2.00) * dxif debug:print "Yn", Ynreturn Ynprint "Result :" , Euler(lambda y: 1+y , 0.0, 1.0 , 1, )