Calculus B -S18.sagews

def lefthand_rs(fcn,a,b,n):
	# output: left-hand riemann sum approx of int_a^n fcn(x)dx using n steps
	Deltax = (b-a)*1.0/n
	return Deltax*sum([f(a+Deltax*i) for i in range(n)])

#Example
n=5
f(x)=e^(-x)
print lefthand_rs(f(x),0,1,n).n();

0.697438279864974

def righthand_rs(fcn,a,b,n):
# output: right-hand riemann sum approx of int_a^b fcn(x)dx using n steps
	Deltax = (b-a)*1.0/n
	return Deltax*sum([f(a+Deltax*(i+1)) for i in range(n)])

#Example
n=5
f(x)=x^2
print righthand_rs(f(x),0,1,n).n();

0.440000000000000

def midpoint_rule(fcn,a,b,n):
# output: midpoint rule approx of int_a^b fcn(x)dx using n steps
	Deltax = (b-a)*1.0/n
	xs=[a+Deltax*i for i in range(n+1)]
	ysmid=[f((xs[i]+xs[i+1])/2) for i in range(n)]
	return Deltax*sum(ysmid)

#Example
n=5;
f(x)=x*x
print midpoint_rule(f(x),0,1,n).n();

0.330000000000000

f(x) = x^2
f = Piecewise([[(0,1), f]])
g = f.riemann_sum(300, mode="right")
F = f.plot(color="blue")
R= add([line([[a,0],[a,f(x=a)],[b,f(x=b)],[b,0]], color="red") for (a,b), f in g.list()])
show(F+R)

def trapezoid_rule(fcn,a,b,n):
# output: trapezoid rule approx of int_a^b fcn(x)dx using n steps
	Deltax = (b-a)*1.0/n
	coeffs = [2]*(n-1)
	coeffs = [1]+coeffs+[1]
	valsf = [f(a+Deltax*i) for i in range(n+1)]
	return (Deltax/2)*sum([coeffs[i]*valsf[i] for i in range(n+1)])

#Example
n=10;
f(x)=1/x
print trapezoid_rule(f(x),1,2,n).n();

0.693771403175428

x = var('x') 
f1(x) = e^(x^2) 
f = Piecewise([[(0,1),f1]]) 
trapezoid_sum = f.trapezoid(5) 
P = f.plot(rgbcolor=(0,0,1), plot_points=40)
Q = trapezoid_sum.plot(rgbcolor=(1,0,0), plot_points=40) 
L = add([line([[a,0],[a,f(x=a)]], rgbcolor=(1,0,0)) for (a,b), f in trapezoid_sum.list()])
M = line([[1,0],[1,f1(1)]], rgbcolor=(1,0,0)) 
show(P + Q + L + M)

def simpsons_rule(fcn,a,b,n):
# output: simpsons rule approx of int_a^b fcn(x)dx using n steps (n must be an even integer)
	Deltax = (b-a)*1.0/n
	n2=int(n/2)
	coeffs = [4,2]*n2
	coeffs = [1] +coeffs[:n-1]+[1]
	valsf = [f(a+Deltax*i) for i in range(n+1)]
	return (Deltax/3)*sum([coeffs[i]*valsf[i] for i in range(n+1)])

#Example
n=7;
f(x)=1/x
print simpsons_rule(f(x),1,2,n).n();

0.668501868501868