SharedCalculus B -S18.sagewsOpen in CoCalc
Sage Codes
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