SharedAM2017.Section.5.1.sagewsOpen in CoCalc
Author: Jiajun Ma
Views : 83
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)
zx=z_x = xye(xy)+e(xy)\displaystyle x y e^{\left(x y\right)} + e^{\left(x y\right)}
zy=z_y = x2e(xy)\displaystyle x^{2} e^{\left(x y\right)}
zxx=z_{xx} = xy2e(xy)+2ye(xy)\displaystyle x y^{2} e^{\left(x y\right)} + 2 \, y e^{\left(x y\right)}
zyy=z_{yy} = x3e(xy)\displaystyle x^{3} e^{\left(x y\right)}
zxy=z_{xy} = x2ye(xy)+2xe(xy)\displaystyle x^{2} y e^{\left(x y\right)} + 2 \, x e^{\left(x y\right)}
zyx=z_{yx} = x2ye(xy)+2xe(xy)\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())
zxx=z_{xx} = 2x2(x2+y2)2+1x2+y2\displaystyle -\frac{2 \, x^{2}}{{\left(x^{2} + y^{2}\right)}^{2}} + \frac{1}{x^{2} + y^{2}}
zyy=z_{yy} = 2y2(x2+y2)2+1x2+y2\displaystyle -\frac{2 \, y^{2}}{{\left(x^{2} + y^{2}\right)}^{2}} + \frac{1}{x^{2} + y^{2}}
zxx+zyy=z_{xx} + z_{yy} = 0\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.041.98=1.04^{1.98}= 1.08075191020342\displaystyle 1.08075191020342
estimate = 1.08000000000000\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())
zxyz_{xy}= xu(x,y)yu(x,y)D0,0(z)(u(x,y),v(x,y))+yu(x,y)xv(x,y)D0,1(z)(u(x,y),v(x,y))+xu(x,y)yv(x,y)D0,1(z)(u(x,y),v(x,y))+xv(x,y)yv(x,y)D1,1(z)(u(x,y),v(x,y))+2xyu(x,y)D0(z)(u(x,y),v(x,y))+2xyv(x,y)D1(z)(u(x,y),v(x,y))\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)) T = Cylindrical('height', ['radius', 'azimuth']) 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) T = Cylindrical('height', ['radius', 'azimuth']) 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) T = Cylindrical('height', ['radius', 'azimuth']) 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])
[[x=0\displaystyle x = 0, y=0\displaystyle y = 0]]
Function value: [0\displaystyle 0]
[[x=2\displaystyle x = \sqrt{2}, y=122\displaystyle y = \frac{1}{2} \, \sqrt{2}, λ=(72)\displaystyle \lambda = \left(-\frac{7}{2}\right)], [x=2\displaystyle x = -\sqrt{2}, y=122\displaystyle y = -\frac{1}{2} \, \sqrt{2}, λ=(72)\displaystyle \lambda = \left(-\frac{7}{2}\right)], [x=1\displaystyle x = 1, y=(1)\displaystyle y = \left(-1\right), λ=1\displaystyle \lambda = 1], [x=(1)\displaystyle x = \left(-1\right), y=1\displaystyle y = 1, λ=1\displaystyle \lambda = 1]]
Function value: [212\displaystyle \frac{21}{2}, 212\displaystyle \frac{21}{2}, 3\displaystyle -3, 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)
23\displaystyle \frac{2}{3} = 23\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))
525\displaystyle \frac{52}{5}
458\displaystyle \frac{45}{8}
sin(1)+1\displaystyle -\sin\left(1\right) + 1
12\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))
32π\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))
19(3π4)R3\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)`
2a2\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))
8936468\displaystyle \frac{893}{6468} = 0.138064316636\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))
0.0142857142857\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)
0.330217428727354\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))
0.9936306176\displaystyle 0.9936306176
0.032773828125\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))
0.00931615460852329\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(X0)=0.000045P(X\leq 0)=0.000045
P(X1)=0.000499P(X\leq 1)=0.000499
P(X2)=0.002769P(X\leq 2)=0.002769
P(X3)=0.010336P(X\leq 3)=0.010336
P(X4)=0.029253P(X\leq 4)=0.029253
P(X5)=0.067086P(X\leq 5)=0.067086
P(X6)=0.130141P(X\leq 6)=0.130141
P(X7)=0.220221P(X\leq 7)=0.220221
P(X8)=0.332820P(X\leq 8)=0.332820
P(X9)=0.457930P(X\leq 9)=0.457930
P(X10)=0.583040P(X\leq 10)=0.583040
P(X11)=0.696776P(X\leq 11)=0.696776
P(X12)=0.791556P(X\leq 12)=0.791556
P(X13)=0.864464P(X\leq 13)=0.864464
P(X14)=0.916542P(X\leq 14)=0.916542
P(X15)=0.951260P(X\leq 15)=0.951260
P(X16)=0.972958P(X\leq 16)=0.972958
P(X17)=0.985722P(X\leq 17)=0.985722
P(X18)=0.992813P(X\leq 18)=0.992813
P(X19)=0.996546P(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))
[μ=0\displaystyle \mu = 0, μ=2\displaystyle \mu = 2]
0.0902235221577418\displaystyle 0.0902235221577418
0.0902235221577\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))
0.0297701813049\displaystyle 0.0297701813049
0.991867757203\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: 0.00306566200976\displaystyle 0.00306566200976
Binomial: 0.00306397596597\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: 0.0369363131138\displaystyle 0.0369363131138
Binomial: 0.0367803266518\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))
512.620626218\displaystyle 512.620626218
0.682689492137\displaystyle 0.682689492137
0.954499736104\displaystyle 0.954499736104
0.997300203937\displaystyle 0.997300203937
0.682689492137\displaystyle 0.682689492137
0.954499736104\displaystyle 0.954499736104
0.997300203937\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 0.9772498680518208\displaystyle 0.9772498680518208 , Route 2 0.9937903346742238\displaystyle 0.9937903346742238
Route 1 0.8413447460685429\displaystyle 0.8413447460685429 , Route 2 0.5\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)) w = plot3d(multpdf1,(x,-3,3),(y,-3,3),adaptive=True, color=rainbow(60, 'rgbtuple'),optical=0.6,frame=False) 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)) w = plot3d(multpdf1,(x,-3,3),(y,-3,3),adaptive=True, color=rainbow(60, 'rgbtuple'),optical=0.6,frame=False) 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)) w = plot3d(multpdf1,(x,-3,3),(y,-3,3),adaptive=True, color=rainbow(100, 'rgbtuple'),optical=0.6,frame=False) 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]))
118\displaystyle \frac{11}{8}
318\displaystyle \frac{31}{8}
74\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)) w = plot3d(30*f(1,1,0),(x,-3,3),(y,-3,3),adaptive=True, color=rainbow(60, 'rgbtuple'),optical=0.6,frame=False) show(xaxis+yaxis+w) xaxis = arrow3d((-4,0,0),(4,0,0)) yaxis = arrow3d((0,-4,0),(0,4,0)) w = plot3d(30*f(1,1,0.5),(x,-3,3),(y,-3,3),adaptive=True, color=rainbow(120, 'rgbtuple'),optical=0.6,frame=False) show(xaxis+yaxis+w) xaxis = arrow3d((-4,0,0),(4,0,0)) yaxis = arrow3d((0,-4,0),(0,4,0)) w = plot3d(30*f(1,1,-0.5),(x,-3,3),(y,-3,3),adaptive=True, color=rainbow(120, 'rgbtuple'),optical=0.6,frame=False) show(xaxis+yaxis+w)
3D rendering not yet implemented
3D rendering not yet implemented
3D rendering not yet implemented
Error in lines 1-1 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 'W' is not defined