Ticket #10983: sol_calculus.py

File sol_calculus.py, 7.7 KB (added by casamayou, 7 years ago)

updated file

Line 
1r"""
2
3sage: var('n k'); p = 4; s = [n + 1]
4(n, k)
5sage: for k in (1..p):
6...     s = s + [factor((((n + 1)^(k + 1) - sum(binomial(k + 1, j)*s[j]
7...              for j in (0..k - 1))) / (k + 1)))]
8sage: s
9[n + 1, 1/2*(n + 1)*n, 1/6*(n + 1)*(2*n + 1)*n, 1/4*(n + 1)^2*n^2, 1/30*(n + 1)*(2*n + 1)*(3*n^2 + 3*n - 1)*n]
10
11"""
12
13
14r"""
15
16sage: var('x h a'); f = function('f')
17(x, h, a)
18sage: g(x) = taylor(f(x), x, a, 3)
19sage: phi(h) = (g(a+3*h)-3*g(a+2*h)+3*g(a+h)-g(a))/h**3
20sage: phi(h).expand()
21D[0, 0, 0](f)(a)
22sage: n = 7; var('x h a'); f = function('f')
23(x, h, a)
24sage: g(x) = taylor(f(x), x, a, n)
25sage: phi(h) = sum(binomial(n,k)*(-1)^(n-k)*g(a+k*h) for k in (0..n))/h**n
26sage: phi(h).expand()
27D[0, 0, 0, 0, 0, 0, 0](f)(a)
28
29"""
30
31
32r"""
33
34sage: theta = 12 * arctan(1/38) + 20 * arctan(1/57) + 7 * arctan(1/239) + 24 * arctan(1/268)
35sage: x = tan(theta)
36sage: y = x.trig_expand()
37sage: y.trig_simplify()
381
39sage: M = 12*(1/38)+20*(1/57)+ 7*(1/239)+24*(1/268)
40sage: M
4137735/48039
42sage: x = var('x')
43sage: f(x) = taylor(arctan(x), x, 0, 21)
44sage: approx = 4 * (12 * f(1/38) + 20 * f(1/57) + 7 * f(1/239) + 24 * f(1/268))
45sage: approx.n(digits = 50); pi.n(digits = 50)
463.1415926535897932384626433832795028851616168852864
473.1415926535897932384626433832795028841971693993751
48sage: approx.n(digits = 50) - pi.n(digits = 50)
499.6444748591132486785420917537404705292978817080880e-37
50
51"""
52
53
54r"""
55
56sage: n = var('n'); phi = lambda x: n*pi +pi/2 - arctan(1/x); x = pi * n
57sage: for i in range(4): x = taylor(phi(x), n, oo, 2 * i); x
581/2*pi + pi*n
591/2*pi + pi*n - 1/(pi*n) + 1/2/(pi*n^2)
601/2*pi + pi*n - 1/(pi*n) + 1/2/(pi*n^2) - 1/12*(3*pi^2 + 8)/(pi^3*n^3) + 1/8*(pi^2 + 8)/(pi^3*n^4)
611/2*pi + pi*n - 1/(pi*n) + 1/2/(pi*n^2) - 1/12*(3*pi^2 + 8)/(pi^3*n^3) + 1/8*(pi^2 + 8)/(pi^3*n^4) - 1/240*(15*pi^4 + 240*pi^2 + 208)/(pi^5*n^5) + 1/96*(3*pi^4 + 80*pi^2 + 208)/(pi^5*n^6)
62
63"""
64
65
66r"""
67
68sage: var('h')
69h
70sage: f(x, y) = x * y * (x**2 - y**2) / (x**2 + y**2)
71sage: D1f(x, y) = diff(f(x,y), x)
72sage: limit((D1f(0,h) - 0) / h, h=0)
73-1
74sage: D2f(x, y) = diff(f(x,y), y)
75sage: limit((D2f(h,0) - 0) / h, h=0)
761
77sage: g = plot3d(f(x, y), (x, -3, 3), (y, -3, 3))
78
79"""
80
81
82r"""
83
84sage: var('n t')
85(n, t)
86sage: v(n) = (4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6))*1/16^n
87sage: assume(8*n+1>0)
88sage: u(n) = integrate((4*sqrt(2)-8*t^3-4*sqrt(2)*t^4-8*t^5)\
89...                     * t^(8*n), t, 0, 1/sqrt(2))
90sage: (u(n)-v(n)).simplify_full()
910
92sage: J = integrate((4*sqrt(2)-8*t^3-4*sqrt(2)*t^4-8*t^5)\
93...                 / (1-t^8), t, 0, 1/sqrt(2))
94sage: J.simplify_full()
95pi + 2*log(sqrt(2) - 1) + 2*log(sqrt(2) + 1)
96sage: ln(exp(J).simplify_log())
97pi
98sage: l = sum(v(n) for n in (0..40)); l.n(digits=60); pi.n(digits=60)
993.14159265358979323846264338327950288419716939937510581474759
1003.14159265358979323846264338327950288419716939937510582097494
101sage: print "%e" % (l-pi).n(digits=60)
102-6.227358e-54
103
104"""
105
106
107r"""
108
109sage: var('X'); ps = lambda f,g : integral(f * g, X, -pi, pi)
110X
111sage: n = 5; Q = sin(X)
112sage: var('a a0 a1 a2 a3 a4 a5'); a= [a0, a1, a2, a3, a4, a5]
113(a, a0, a1, a2, a3, a4, a5)
114sage: P = sum(a[k] * X^k for k in (0..n))
115sage: equ = [ps(P - Q, X^k) for k in (0..n)]
116sage: sol = solve(equ, a)
117sage: P = sum(sol[0][k].rhs() * X^k for k in (0..n))
118sage: g = plot(P,X,-4,4,color='red') + plot(Q,X,-4,4,color='blue')
119
120"""
121
122
123r"""
124
125sage: var('p e theta1 theta2 theta3')
126(p, e, theta1, theta2, theta3)
127sage: r(theta) = p / (1-e * cos(theta))
128sage: r1 = r(theta1); r2 = r(theta2); r3 = r(theta3)
129sage: R1 = vector([r1 * cos(theta1), r1 * sin(theta1), 0])
130sage: R2 = vector([r2 * cos(theta2), r2 * sin(theta2), 0])
131sage: R3 = vector([r3 * cos(theta3), r3 * sin(theta3), 0])
132sage: D = R1.cross_product(R2) + R2.cross_product(R3) + R3.cross_product(R1)
133sage: i = vector([1, 0, 0])
134sage: S = (r1 - r3) * R2 + (r3 - r2) * R1 +   (r2 - r1) * R3
135sage: V =  S + e * i.cross_product(D)
136sage: map(lambda x:x.simplify_full(), V) # rep. : [0, 0, 0]
137[0, 0, 0]
138sage: map(lambda x:x.simplify_full(), S.cross_product(D))
139[(e*p^4*sin(theta1)^2*cos(theta2)^2 - 2*e*p^4*sin(theta1)*sin(theta2)*cos(theta1)*cos(theta2) + e*p^4*sin(theta2)^2*cos(theta1)^2 + (e*p^4*cos(theta1)^2 - 2*e*p^4*cos(theta1)*cos(theta2) + e*p^4*cos(theta2)^2)*sin(theta3)^2 + (e*p^4*sin(theta1)^2 - 2*e*p^4*sin(theta1)*sin(theta2) + e*p^4*sin(theta2)^2)*cos(theta3)^2 - 2*(e*p^4*sin(theta1)^2*cos(theta2) + e*p^4*sin(theta2)^2*cos(theta1) - (e*p^4*sin(theta1)*cos(theta1) + e*p^4*sin(theta1)*cos(theta2))*sin(theta2))*cos(theta3) + 2*(e*p^4*sin(theta1)*cos(theta1)*cos(theta2) - e*p^4*sin(theta1)*cos(theta2)^2 - (e*p^4*cos(theta1)^2 - e*p^4*cos(theta1)*cos(theta2))*sin(theta2) - (e*p^4*sin(theta1)*cos(theta1) - e*p^4*sin(theta1)*cos(theta2) - (e*p^4*cos(theta1) - e*p^4*cos(theta2))*sin(theta2))*cos(theta3))*sin(theta3))/(e^2*cos(theta1)^2 + (e^4*cos(theta1)^2 - 2*e^3*cos(theta1) + e^2)*cos(theta2)^2 + (e^4*cos(theta1)^2 - 2*e^3*cos(theta1) + (e^6*cos(theta1)^2 - 2*e^5*cos(theta1) + e^4)*cos(theta2)^2 - 2*(e^5*cos(theta1)^2 - 2*e^4*cos(theta1) + e^3)*cos(theta2) + e^2)*cos(theta3)^2 - 2*(e^3*cos(theta1)^2 - 2*e^2*cos(theta1) + e)*cos(theta2) - 2*(e^3*cos(theta1)^2 + (e^5*cos(theta1)^2 - 2*e^4*cos(theta1) + e^3)*cos(theta2)^2 - 2*e^2*cos(theta1) - 2*(e^4*cos(theta1)^2 - 2*e^3*cos(theta1) + e^2)*cos(theta2) + e)*cos(theta3) - 2*e*cos(theta1) + 1), 0, 0]
140sage: N = r3 * R1.cross_product(R2) + r1 * R2.cross_product(R3) + r2 * R3.cross_product(R1)
141sage: W =  p * S + e * i.cross_product(N)
142sage: print map(lambda x:x.simplify_full(), W)  # rep. : [0, 0, 0]
143[0, 0, 0]
144sage: R1=vector([0,1.,0]);R2=vector([2.,2.,0]);R3=vector([3.5,0,0])
145sage: r1 = R1.norm(); r2 = R2.norm(); r3 = R3.norm()
146sage: D = R1.cross_product(R2) + R2.cross_product(R3) + R3.cross_product(R1)
147sage: S = (r1 - r3) * R2 + (r3 - r2) * R1 + (r2 - r1) * R3
148sage: V =  S + e * i.cross_product(D)
149sage: N = r3 * R1.cross_product(R2) + r1 * R2.cross_product(R3) \
150...  + r2 * R3.cross_product(R1)
151sage: W =  p * S + e * i.cross_product(N)
152sage: e = S.norm() / D.norm()
153sage: p = N.norm() / D.norm()
154sage: a = p/(1-e^2)
155sage: c = a * e
156sage: b = sqrt(a^2 - c^2)
157sage: X = S.cross_product(D)
158sage: i = X / X.norm()
159sage: phi = atan2(i[1],i[0]) * 180 / pi.n()
160sage: print "%.3f %.3f %.3f %.3f %.3f %.3f" % (a, b, c, e, p, phi)
1612.360 1.326 1.952 0.827 0.746 17.917
162
163"""
164
165
166r"""
167
168sage: A = matrix(QQ, [[2, -3, 2, -12, 33],
169...                   [ 6, 1, 26, -16, 69],
170...                   [10, -29, -18, -53, 32],
171...                   [2, 0, 8, -18, 84]])
172sage: A.right_kernel()
173Vector space of degree 5 and dimension 2 over Rational Field
174Basis matrix:
175[     1      0  -7/34   5/17   1/17]
176[     0      1  -3/34 -10/17  -2/17]
177sage: H = A.echelon_form()
178sage: A.column_space()
179Vector space of degree 4 and dimension 3 over Rational Field
180Basis matrix:
181[       1        0        0 1139/350]
182[       0        1        0    -9/50]
183[       0        0        1   -12/35]
184sage: S.<x,y,z,t>=QQ[]
185sage: C = matrix(S, 4,1,[x,y,z,t])
186sage: B = block_matrix([A,C], ncols=2)
187sage: C = B.echelon_form()
188sage: C[3,5]*350
189-1139*x + 63*y + 120*z + 350*t
190sage: K = A.kernel(); K
191Vector space of degree 4 and dimension 1 over Rational Field
192Basis matrix:
193[        1  -63/1139 -120/1139 -350/1139]
194sage: matrix(K.0).right_kernel()
195Vector space of degree 4 and dimension 3 over Rational Field
196Basis matrix:
197[       1        0        0 1139/350]
198[       0        1        0    -9/50]
199[       0        0        1   -12/35]
200sage: A = matrix(QQ, [[-2, 1, 1], [8, 1, -5], [4, 3, -3]])
201sage: C = matrix(QQ, [[1, 2, -1], [2, -1, -1], [-5, 0, 3]])
202sage: B = C.solve_left(A); B
203[ 0 -1  0]
204[ 2  3  0]
205[ 2  1  0]
206sage: C.left_kernel()
207Vector space of degree 3 and dimension 1 over Rational Field
208Basis matrix:
209[1 2 1]
210sage: var('x y z'); v = matrix([[1, 2, 1]])
211(x, y, z)
212sage: B = B+(x*v).stack(y*v).stack(z*v); B
213[      x 2*x - 1       x]
214[  y + 2 2*y + 3       y]
215[  z + 2 2*z + 1       z]
216sage: A == B*C
217True
218
219
220"""
221