Shared2017-09-19-120509.sagewsOpen in CoCalc
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=10;
f(x)=1/x
print simpsons_rule(f(x),1,2,n).n();
0.693150230688930



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



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=10;
f(x)=1/x
print midpoint_rule(f(x),1,2,n).n();
0.692835360409960


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=10
f(x)=1/x
print righthand_rs(f(x),1,2,n).n();
0.668771403175428

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)/n
	return Deltax*sum([f(a+Deltax*i) for i in range(n)])

#Example
n=6
f(x)=1/x
print lefthand_rs(f(x),1,2,n).n();
0.736544011544012

# the actual integral gives
f(x)=1/x
integral(f , (x,1,2))
N(log(2))
log(2) 0.693147180559945



x = var('x') 
f1(x) = 1/x
f = Piecewise([[(1,2),f1]]) 
trapezoid_sum = f.trapezoid(10) 
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([[2,0],[2,f1(2)]], rgbcolor=(1,0,0)) 
show(P + Q + L + M)


f(x) = x^2-5*x+10
f = Piecewise([[(0,10), f]])
g = f.riemann_sum(6, mode="midpoint")
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)

integral(1/(x^2+1),(x, 0, oo)) # we getting pi/2
plot(1/(x^2+1),x,0,200) # this gives us that the idea that integral of the function will converge from 0 to infinity.
1/2*pi

plot(1/(sqrt(x-2)),2,5) # discontinuous integrands (Improper integral of type 2)
integral(1/(sqrt(x-2)),2,5) # this integral is convergent

2*sqrt(3)


f(x)=1/(sqrt(x^2-x))
g(x)=1/x
plot(f, 3, 7,gridlines='minor', frame=True, color='red')+ plot(g,3 ,7)


f(x)=x/(sqrt(x^2+1))
g(x)=x^4-x
plot(f, 0, 1.4,gridlines='minor', frame=True)+ plot(g,0 ,1.4)+ text('$\int_0^{1.18}(f-g)(x)\,dx$',(0.5,2),fontsize=20)



f(x)=x^2
g(x)=2*x-x^2
plot(f(x), 0, 2)+ plot(g(x),0 ,2)
integral(g(x)-f(x),(x,0,1))
q(x)=f(x)-g(x)
solve(q(x),x)
1/3 [x == 0, x == 1]

f(x)=e^(x)
g(x)=x
plot(f(x), 0, 1,color='green', legend_label='f(x)=e^(x)')+ plot(g(x),0 ,1,color='red', legend_label='g(x)=x')
integral(f(x)-g(x),(x,0,1))
e - 3/2


f(x)=x/(sqrt(x^2+1))
g(x)=x^4-x
plot(f(x), 0, 1.4,gridlines=True, frame=False)+ plot(g(x),0 ,1.4)
integral(f(x)-g(x),(x,0,1.18))
q(x)=f(x)-g(x)
solve(q(x),x)
find_root(q(x), 1, 1.2)
0.7853869527188246 [x == 0, x^3 == (sqrt(x^2 + 1) + 1)/sqrt(x^2 + 1)] 1.1807757031062647

solve(x^3-7*x^2+14*x-8, x)
[x == 2, x == 4, x == 1]


find_root(cos(x)-x^2+1, 1,2)
1.1765019399018324

f(x)=cos(x)
g(x)=x^2-1
plot(f, -2, 2,gridlines='minor', frame=True, color='red')+ plot(g, -2 ,2)
integral(2*(f-g),0,1.18)
3.113857358149374

plot(atan(x)/(2+e^x),x)

s = 'sum (5*(-2/3)^(n-1),n,1,inf), simpsum'
SR(sage.calculus.calculus.maxima(s))
3

s = 'sum (1/n^2,n,1,inf), simpsum'
SR(sage.calculus.calculus.maxima(s))
1/6*pi^2

s = 'sum (1/2^n,n,1,inf), simpsum'
SR(sage.calculus.calculus.maxima(s))
1


s = 'sum (1/2^n,n,1,5), simpsum'
SR(sage.calculus.calculus.maxima(s))
numerical_approx(maxima(s), digits=10)
31/32 0.9687500000




print " Plot Slope Field "
t , v = var ("t,v")
pt = plot_slope_field ( 9.8-0.2*v , (t , 0 , 10) ,(v , 40 , 55) , plot_points = \
20)
show ( pt , aspect_ratio ='0.5')
Plot Slope Field

t = var('t')
x = function('x', t)
de = lambda v: diff(v,t) == 9.8-0.2*v
desolve(de(x),[x,t]);
(_C + 49*e^(1/5*t))*e^(-1/5*t)

print " Plot Slope Field "
t , p = var ("t,p")
pt = plot_slope_field ( 0.5*p-450 , (t , 0 , 10) ,(p , 800 , 1000) , plot_points = \
20)
show ( pt , aspect_ratio ='0.02')
Plot Slope Field


t = var('t')
x = function('x', t)
de = lambda p: diff(p,t) == 0.5*p-450
desolve(de(x),[x,t]);

integrate(e^x/(3-e^(2*x)),x)
-1/6*sqrt(3)*log(-(sqrt(3) - e^x)/(sqrt(3) + e^x))

plot(e^x/(3-e^(2*x)),(x,0.54,0.56))

integrate(1/(x**2),(x,1,oo))
1