CoCalc Public FilesAM2017.Section.5.1.sagews
Author: Jiajun Ma
Views : 205
Description: 高等数学（交大生农医药2017版）Section 5.1
# Example 5.1.5
x, y,t,r = var('x y t r')
def cf(r,t): return 2*sin(t)*cos(t)*r
P = parametric_plot3d((r*sin(t),r*cos(t),sin(t)*cos(t)),(r,0,.7),(t,-pi,pi),color=(cf,colormaps.ocean),opacity=0.7,frame=false)
PL=plot3d(0, (-1,1), (-1,1),opacity=0.2,mash=false,frame=false)
from sage.plot.plot3d.plot3d import axes
T = axes(1, .5);
(P+PL+T).show()

3D rendering not yet implemented
# Example 5.1.10
x,y = var('x y')
z = x*exp(x*y);
z_x = diff(z,x);show('$z_x =$',z_x)
z_y = diff(z,y);show('$z_y =$',z_y)
z_xx = diff(z_x,x);show('$z_{xx} =$',z_xx)
z_yy = diff(z_y,y);show('$z_{yy} =$',z_yy)
z_xy = diff(z_y,x);show('$z_{xy} =$',z_xy)
z_yx = diff(z_x,y);show('$z_{yx} =$',z_yx)


$z_x =$ $\displaystyle x y e^{\left(x y\right)} + e^{\left(x y\right)}$
$z_y =$ $\displaystyle x^{2} e^{\left(x y\right)}$
$z_{xx} =$ $\displaystyle x y^{2} e^{\left(x y\right)} + 2 \, y e^{\left(x y\right)}$
$z_{yy} =$ $\displaystyle x^{3} e^{\left(x y\right)}$
$z_{xy} =$ $\displaystyle x^{2} y e^{\left(x y\right)} + 2 \, x e^{\left(x y\right)}$
$z_{yx} =$ $\displaystyle x^{2} y e^{\left(x y\right)} + 2 \, x e^{\left(x y\right)}$
# Example 5.1.11
x,y = var('x y')
z = ln(sqrt(x^2+y^2));
z_xx = diff(z,x,2);show('$z_{xx} =$',z_xx)
z_yy = diff(z,y,2);show('$z_{yy} =$',z_yy)
show('$z_{xx} + z_{yy} =$', (z_xx+z_yy).simplify_full())

$z_{xx} =$ $\displaystyle -\frac{2 \, x^{2}}{{\left(x^{2} + y^{2}\right)}^{2}} + \frac{1}{x^{2} + y^{2}}$
$z_{yy} =$ $\displaystyle -\frac{2 \, y^{2}}{{\left(x^{2} + y^{2}\right)}^{2}} + \frac{1}{x^{2} + y^{2}}$
$z_{xx} + z_{yy} =$ $\displaystyle 0$
# Example estimate 1.04^{1.98}
show("$1.04^{1.98}=$",exp(ln(1.04)*1.98))
z(x,y) = x^y;
x0=1;y0=2;
dx = 0.04; dy=-0.02;
show("estimate = ", z(x0,y0)+ diff(z,x)(x0,y0)*dx+diff(z,y)(x0,y0)*dy)

$1.04^{1.98}=$ $\displaystyle 1.08075191020342$
estimate = $\displaystyle 1.08000000000000$
# Section 5.1.3.3
x, y = var('x y')
u=function('u')(x,y)
v=function('v')(x,y)
z=function('z')(u,v)
show("$z_{xy}$=",diff(diff(z,x),y).expand())

$z_{xy}$= $\displaystyle \frac{\partial}{\partial x}u\left(x, y\right) \frac{\partial}{\partial y}u\left(x, y\right) \mathrm{D}_{0, 0}\left(z\right)\left(u\left(x, y\right), v\left(x, y\right)\right) + \frac{\partial}{\partial y}u\left(x, y\right) \frac{\partial}{\partial x}v\left(x, y\right) \mathrm{D}_{0, 1}\left(z\right)\left(u\left(x, y\right), v\left(x, y\right)\right) + \frac{\partial}{\partial x}u\left(x, y\right) \frac{\partial}{\partial y}v\left(x, y\right) \mathrm{D}_{0, 1}\left(z\right)\left(u\left(x, y\right), v\left(x, y\right)\right) + \frac{\partial}{\partial x}v\left(x, y\right) \frac{\partial}{\partial y}v\left(x, y\right) \mathrm{D}_{1, 1}\left(z\right)\left(u\left(x, y\right), v\left(x, y\right)\right) + \frac{\partial^{2}}{\partial x\partial y}u\left(x, y\right) \mathrm{D}_{0}\left(z\right)\left(u\left(x, y\right), v\left(x, y\right)\right) + \frac{\partial^{2}}{\partial x\partial y}v\left(x, y\right) \mathrm{D}_{1}\left(z\right)\left(u\left(x, y\right), v\left(x, y\right)\right)$
# Example 5.1.5
x, y = var('x y')
def cf(x,y): return abs(x*y)
P = plot3d(x*y,(x,-1,1),(y,-1,1),color=(cf,colormaps.rainbow),opacity=0.7,frame=false)
from sage.plot.plot3d.plot3d import axes
T = axes(1, .5,color='black');
(P+T).show()



# Example 5.1.5-1
def cf(r,t): return abs(3*r^2*cos(t)*sin(t))
r, t, z = var('r theta z')
#P=plot3d(r^2, (r, 0, 1), (t, -pi, pi), color=(cf,colormaps.rainbow),opacity=0.7,frame=false)
PP = plot3d(r^2*cos(t)*sin(t), (r, 0, 1), (t, 0, 2*pi), transformation=T,color=(cf,colormaps.rainbow),axes=true, adaptive=True,frame=false,opacity=0.7,aspect_ratio=[1,1,1.5])
from sage.plot.plot3d.plot3d import axes
CD = axes(1, .5,color='black');
(PP+CD).show()

3D rendering not yet implemented
def cf(r,theta): return abs(r^2)
r, t, z = var('r theta z')
#P=plot3d(r^2, (r, 0, 1), (t, -pi, pi), color=(cf,colormaps.rainbow),opacity=0.7,frame=false)
PP = plot3d(r^2, (r, 0, 1), (t, 0, 2*pi), transformation=T,color=(cf,colormaps.rainbow),axes=true, frame=false,opacity=0.7,aspect_ratio=[1,1,1.5])
from sage.plot.plot3d.plot3d import axes
CD = axes(1, .5,color='black');
(PP+CD).show()

def cf(r,theta): return abs(r^2)
r, t, z = var('r theta z')
#P=plot3d(r^2, (r, 0, 1), (t, -pi, pi), color=(cf,colormaps.rainbow),opacity=0.7,frame=false)
PP = plot3d(2*r^2, (r, 0, 0.5), (t, 0, 2*pi), transformation=T,color=(cf,colormaps.rainbow),axes=true,
frame=false,opacity=0.7,aspect_ratio=1)

def cff(r,t): return cos(4*r)
PPP = plot3d((1-cos(4*r))/4, (r, 0, 1), (t, 0, 1.5*pi), transformation=T,
#color=(cff,colormaps.rainbow),mesh=true,
frame=false,opacity=0.7,aspect_ratio=1)
from sage.plot.plot3d.plot3d import axes
CD = axes(1.5, .5,color='black');
(PP+PPP+CD).show()

3D rendering not yet implemented
# Example 5.1.30
x,y,lam= var('x y lambda_');
F = 2*x^2+6*x*y+y^2
S0=solve([diff(F,x)==0,diff(F,y)==0],x,y)
show(S0)
show("Function value:",[F.subs(s) for s in S0])

G = x^2+2*y^2-3
L = F + lam*G
Lx = diff(L,x)
Ly = diff(L,y)
Llam = diff(L,lam)
S = solve([Lx==0,Ly==0,Llam==0],x,y,lam)
show(S)
show("Function value:",[F.subs(s) for s in S])

[[$\displaystyle x = 0$, $\displaystyle y = 0$]]
Function value: [$\displaystyle 0$]
[[$\displaystyle x = \sqrt{2}$, $\displaystyle y = \frac{1}{2} \, \sqrt{2}$, $\displaystyle \lambda = \left(-\frac{7}{2}\right)$], [$\displaystyle x = -\sqrt{2}$, $\displaystyle y = -\frac{1}{2} \, \sqrt{2}$, $\displaystyle \lambda = \left(-\frac{7}{2}\right)$], [$\displaystyle x = 1$, $\displaystyle y = \left(-1\right)$, $\displaystyle \lambda = 1$], [$\displaystyle x = \left(-1\right)$, $\displaystyle y = 1$, $\displaystyle \lambda = 1$]]
Function value: [$\displaystyle \frac{21}{2}$, $\displaystyle \frac{21}{2}$, $\displaystyle -3$, $\displaystyle -3$]
# Example on ppt page 17
x,y = var('x y')
I1 = integral(integral(y,y,0,sqrt(1-x^2)),x,-1,1)
I2 = integral(integral(y+y^2*sin(x),y,0,sqrt(1-x^2)),x,-1,1)
show(I1,"=",I2)

$\displaystyle \frac{2}{3}$ = $\displaystyle \frac{2}{3}$
# Example on ppt page 15
x, y = var('x y')
show(integral(integral(x+2*y,y,0,x^2),x,0,2))
show(integral(integral(x*y,x,y^2,y+2),y,-1,2))

assume(x>0);assume(y>0);
# show(integral(integral(sin(y)/y,y,x,sqrt(x)),x,0,1)) # cause error!
show(integral(integral(sin(y)/y,x,y^2,y),0,1))

#Example on ppt page 16
# show(integral(integral(cos(y)/y,y,x,pi/6),x,0,pi/6)) # cause error!
show(integral(integral(cos(y)/y,x,0,y),y,0,pi/6))

$\displaystyle \frac{52}{5}$
$\displaystyle \frac{45}{8}$
$\displaystyle -\sin\left(1\right) + 1$
$\displaystyle \frac{1}{2}$
# Example on ppt page 21
r,t=var('r t')
show(integral(integral(r^3,r,0,2*cos(t)),t,-pi/2,pi/2))

$\displaystyle \frac{3}{2} \, \pi$
# Example on ppt page 22
r,t,R=var('r t R')
assume(R>0)
#show(integral(r*sqrt(R^2-r^2),r,0,R*cos(t)).expand())
show(2*integral(1/3*R^3*(1-(sin(t))^3),t,0,pi/2))

$\displaystyle \frac{1}{9} \, {\left(3 \, \pi - 4\right)} R^{3}$
# Example on ppt page 23
r, t,a = var('r t a')
assume(a>0)
A = 4*integral(integral(r,r,0,sqrt(2*a^2*((cos(t))^2-(sin(t))^2))),t,0,pi/4)
show(A)

$\displaystyle 2 \, a^{2}$
#################
### Chapter 6 ###
#################




# Example 6.2.1
P = binomial(100-3,4)*binomial(3,1)/binomial(100,5);
show(P,'=',float(P))

$\displaystyle \frac{893}{6468}$ = $\displaystyle 0.138064316636$
# Example
P = binomial(50-3,2)*binomial(3,2)/binomial(50,4)+binomial(50-3,1)*binomial(3,3)/binomial(50,4)
show(float(P))

$\displaystyle 0.0142857142857$
# Example 6.2.3
def Pc(n):
return factorial(n)*binomial(365,n)/365^n
for i in xrange(2,51):
print i, 'students; ', float(1-Pc(i))


2 students; 0.0027397260274 3 students; 0.00820416588478 4 students; 0.0163559124666 5 students; 0.0271355736998 6 students; 0.0404624836491 7 students; 0.056235703096 8 students; 0.0743352923517 9 students; 0.0946238338892 10 students; 0.116948177711 11 students; 0.141141378322 12 students; 0.167024788838 13 students; 0.194410275232 14 students; 0.223102512005 15 students; 0.252901319764 16 students; 0.283604005253 17 students; 0.315007665297 18 students; 0.346911417872 19 students; 0.379118526032 20 students; 0.411438383581 21 students; 0.443688335165 22 students; 0.475695307663 23 students; 0.507297234324 24 students; 0.538344257915 25 students; 0.568699703969 26 students; 0.598240820136 27 students; 0.626859282263 28 students; 0.654461472342 29 students; 0.680968537478 30 students; 0.706316242719 31 students; 0.730454633729 32 students; 0.75334752785 33 students; 0.774971854176 34 students; 0.79531686462 35 students; 0.814383238875 36 students; 0.83218210638 37 students; 0.848734008216 38 students; 0.864067821082 39 students; 0.878219664367 40 students; 0.891231809818 41 students; 0.903151611482 42 students; 0.914030471562 43 students; 0.923922855656 44 students; 0.932885368551 45 students; 0.940975899466 46 students; 0.948252843367 47 students; 0.954774402833 48 students; 0.960597972879 49 students; 0.965779609323 50 students; 0.970373579578
# Example on ppt page 40
show(1-(1-0.004)^100)

$\displaystyle 0.330217428727354$
# Example on ppt page 43
def b(k,n,p):
return binomial(n,k)*p^k*(1-p)^(n-k)
P = sum([b(i,10,12/60) for i in xrange(0,6)])
show(float(P))

# Example 6.3.8
P = 1-b(0,6,1/20)-b(1,6,1/20);
show(float(P))

$\displaystyle 0.9936306176$
$\displaystyle 0.032773828125$
# Example on ppt page 47
PA = 0.0004
PB_A=0.94
PB_AC=0.04
show(PA*PB_A/(PA*PB_A+(1-PA)*PB_AC))

$\displaystyle 0.00931615460852329$
# Example on 6.5 ppt page 8
from scipy.stats import binom
binom_dist = binom(5,0.1)
# bar_chart([binom_dist.pmf(x) for x in range(6)])
for x in range(6):
print "P(X=%d) = %f"%(x,binom_dist.pmf(x))

P(X=0) = 0.590490 P(X=1) = 0.328050 P(X=2) = 0.072900 P(X=3) = 0.008100 P(X=4) = 0.000450 P(X=5) = 0.000010
# Example on 6.5 ppt page 9
from scipy.stats import binom
binom_dist = binom(8,0.3)
# bar_chart([binom_dist.pmf(x) for x in range(6)])
for x in range(9):
print "P(X=%d) = %f"%(x,binom_dist.pmf(x))

print "P(X>=2) = %f"%(1-binom_dist.pmf(0)-binom_dist.pmf(1))   # or use: (1-binom_dist.cdf(1))

P(X=0) = 0.057648 P(X=1) = 0.197650 P(X=2) = 0.296475 P(X=3) = 0.254122 P(X=4) = 0.136137 P(X=5) = 0.046675 P(X=6) = 0.010002 P(X=7) = 0.001225 P(X=8) = 0.000066 P(X>=2) = 0.744702
# Example 6.5.3
from scipy.stats import poisson
for k in xrange(0,20):
show("$P(X\leq %i)=%f$"%(k, poisson.cdf(k, mu=10)))

$P(X\leq 0)=0.000045$
$P(X\leq 1)=0.000499$
$P(X\leq 2)=0.002769$
$P(X\leq 3)=0.010336$
$P(X\leq 4)=0.029253$
$P(X\leq 5)=0.067086$
$P(X\leq 6)=0.130141$
$P(X\leq 7)=0.220221$
$P(X\leq 8)=0.332820$
$P(X\leq 9)=0.457930$
$P(X\leq 10)=0.583040$
$P(X\leq 11)=0.696776$
$P(X\leq 12)=0.791556$
$P(X\leq 13)=0.864464$
$P(X\leq 14)=0.916542$
$P(X\leq 15)=0.951260$
$P(X\leq 16)=0.972958$
$P(X\leq 17)=0.985722$
$P(X\leq 18)=0.992813$
$P(X\leq 19)=0.996546$
# Example 2 on ppt page 13
from math import factorial,exp
from scipy.stats import poisson
mu=var("mu")
S = solve([mu/1==mu^2/2],mu)
show(S)
show(2.0^4/(factorial(4))*exp(-2.0))
show(poisson.pmf(k=4,mu=2))

[$\displaystyle \mu = 0$, $\displaystyle \mu = 2$]
$\displaystyle 0.0902235221577418$
$\displaystyle 0.0902235221577$
# Example 3 on ppt page 13
from scipy.stats import poisson
show(poisson.pmf(k=8, mu=4))
show(poisson.cdf(k=9,mu=4))

$\displaystyle 0.0297701813049$
$\displaystyle 0.991867757203$
# Example on ppt page 14-2
from scipy.stats import poisson,binom
show("Poisson: \t\t",poisson.pmf(5, mu=1))
show("Binomial:\t\t", binom.pmf(5,10000,0.0001))

Poisson: $\displaystyle 0.00306566200976$
Binomial: $\displaystyle 0.00306397596597$
# Example on ppt page 14-3
from scipy.stats import poisson,binom
show("Poisson:",1-poisson.cdf(1, mu=150.0/500))
show("Binomial:", 1-binom.cdf(1,150,1.0/500))

Poisson: $\displaystyle 0.0369363131138$
Binomial: $\displaystyle 0.0367803266518$
from scipy.stats import binom
from sage.plot.point import point
from sage.plot.line import line
n=5
p=0.1
t = text("cdf of B(%i,%g)"%(n,p), (n,1.1), rgbcolor=(1,0,0))
P = [(x,binom.cdf(x,n,p)) for x in range(-1,n+1)]
t = t+ point(P[1:],size=50)
for x, y in P:
t = t+ line([(x,y),(x+1,y)],thickness=2)
t.show()

from scipy.stats import norm
p=plot(norm.cdf,xmin=-3.5,xmax=3.5)
show(p+text("$\Phi$", (0.2,1)))
q=plot(norm.pdf,xmin=-3.5,xmax=3.5,rgbcolor=(0,.5,1))
show(q)
show(p+q)

# Example on Section 6.5 ppt p33
from scipy.stats import norm
show(400*norm.ppf(0.9))

show(1-2*norm.cdf(-1))
show(1-2*norm.cdf(-2))
show(1-2*norm.cdf(-3))


$\displaystyle 512.620626218$
$\displaystyle 0.682689492137$
$\displaystyle 0.954499736104$
$\displaystyle 0.997300203937$
$\displaystyle 0.682689492137$
$\displaystyle 0.954499736104$
$\displaystyle 0.997300203937$
# Example on Section 6.5 ppt 6
from scipy.stats import norm
show("Route 1",norm.cdf(70,loc=50,scale=10),", Route 2",norm.cdf(70,loc=60,scale=4))
show("Route 1",norm.cdf(60,loc=50,scale=10),", Route 2",norm.cdf(60,loc=60,scale=4))


Route 1 $\displaystyle 0.9772498680518208$ , Route 2 $\displaystyle 0.9937903346742238$
Route 1 $\displaystyle 0.8413447460685429$ , Route 2 $\displaystyle 0.5$
# 2d normal distribution
from sage.plot.contour_plot import ContourPlot

x,y,u1, u2, s1,s2,rho=var("x y u1 u2 s1 s2 rho")
#def mvpdf(x,y,u1,u2,s1,s2,rho):
#    assert(abs(rho)<1)
#   return
multpdf = 1/(2*pi*sqrt(1-rho^2))*exp(-0.5/(1-rho^2)*((x-u1)^2/s1^2+(y-u2)^2/s2^2-2*rho*(x-u1)*(y-u2)/s1/s2))

multpdf1 = 30*multcdf.subs(u1=0,u2=0,s1=1,s2=1,rho=0)
xaxis = arrow3d((-4,0,0),(4,0,0))
yaxis = arrow3d((0,-4,0),(0,4,0))
show(xaxis+yaxis+w)

multpdf1 = 30*multcdf.subs(u1=0,u2=0,s1=1,s2=1,rho=0.8)
xaxis = arrow3d((-4,0,0),(4,0,0))
yaxis = arrow3d((0,-4,0),(0,4,0))
show(xaxis+yaxis+w)

multpdf1 = 30*multcdf.subs(u1=0,u2=0,s1=1,s2=1,rho=-0.8)
xaxis = arrow3d((-4,0,0),(4,0,0))
yaxis = arrow3d((0,-4,0),(0,4,0))
show(xaxis+yaxis+w)

Error in lines 3-3 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute flags=compile_flags) in namespace, locals File "", line 1, in <module> NameError: name 'multcdf' is not defined

k = var('k')
g = plot(exp(-ln(k)/k),(k,2,30))
show(g)


show((2*120+2*130+3*125+110+135+140.0)/10)
D1 = [0.028, 0.121,0.234,0.267,0.200,0.103,0.037,0.009,0.001]
D2 = [0.001, 0.010,0.044, 0.117, 0.205, 0.247, 0.205, 0.117, 0.044, 0.010]
show(sum([i*p for i,p in enumerate(D1)]))
show(sum([i*p for i,p in enumerate(D2)]))

D = [(-1,1/8),(0,1/4),(2,3/8),(3,1/4)]
show(sum([x*p for x,p in D]))
show(sum([x^2 *p for x,p in D]))
show(sum([(-2*x+1) *p for x,p in D]))

$\displaystyle \frac{11}{8}$
$\displaystyle \frac{31}{8}$
$\displaystyle -\frac{7}{4}$



x, y = var('x y')
def f(s,t,r):
return 1/(2*pi*s*t*sqrt(1-r^2))*exp(-(x^2/(s^2)-2*r*x*y/s/t+y^2/(t^2))/(2*(1-r^2)))

xaxis = arrow3d((-4,0,0),(4,0,0))
yaxis = arrow3d((0,-4,0),(0,4,0))
show(xaxis+yaxis+w)

xaxis = arrow3d((-4,0,0),(4,0,0))
yaxis = arrow3d((0,-4,0),(0,4,0))
show(xaxis+yaxis+w)

xaxis = arrow3d((-4,0,0),(4,0,0))
yaxis = arrow3d((0,-4,0),(0,4,0))
show(xaxis+yaxis+w)

3D rendering not yet implemented
3D rendering not yet implemented
3D rendering not yet implemented
from scipy.stats import binom
from sage.plot.point import point
from sage.plot.line import line
n=2
p=0.6
t = text("cdf of B(%i,%g)"%(n,p), (n,1.1), rgbcolor=(1,0,0))
P = [(x,binom.cdf(x,n,p)) for x in range(-1,n+1)]
for x, y in P:
t = t+ line([(x,y),(x+1,y)],thickness=2,rgbcolor=(0,0,1))
t = t+ point(P[1:],size=50)
for x,y in P[:-1]:
t=t+circle((x+1,y),0.03,fill=True,facecolor=(1,1,1))
t.show()



`