Ticket #10983: calculus.py

File calculus.py, 10.1 KB (added by casamayou, 6 years ago)

updated file

Line 
1r"""
2
3sage: var ('a, x') ; y = cos(x+a)*(x+1) ; y
4(a, x)
5(x + 1)*cos(a + x)
6sage: var ('a,x') ; y = cos(x+a)*(x+1) ; y
7(a, x)
8(x + 1)*cos(a + x)
9sage: y.subs(a=-x) ; y.subs(x=pi/2,a=pi/3) ; y.subs(x=0.5,a=2.3)
10x + 1
11-1/4*(pi + 2)*sqrt(3)
12-1.41333351100299
13sage: y(a=-x) ; y(x=pi/2,a=pi/3) ; y(x=0.5,a=2.3)
14x + 1
15-1/4*(pi + 2)*sqrt(3)
16-1.41333351100299
17sage: var('x,y,z') ; q = x*y+y*z+z*x
18(x, y, z)
19sage: bool (q(x=y,y=z,z=x)==q), bool(q(z=y)(y=x) == 3*x^2)
20(True, True)
21
22"""
23
24
25r"""
26
27sage: var('y z'); f = x^3+y^2+z; f.subs_expr(x^3==y^2, z==1)
28(y, z)
292*y^2 + 1
30sage: f(x)=(2*x+1)^3 ; f(-3)
31  -125
32sage: f(x).expand()
338*x^3 + 12*x^2 + 6*x + 1
34sage: var('y'); u = sin(x) + x*cos(y)
35y
36sage: v = u.function(x, y); v
37(x, y) |--> x*cos(y) + sin(x)
38sage: w(x, y) = u; w
39(x, y) |--> x*cos(y) + sin(x)
40
41"""
42
43
44r"""
45
46sage: (x^x/x).simplify()
47x^(x - 1)
48sage: f = (e^x-1)/(1+e^(x/2)); f.simplify_exp()
49e^(1/2*x) - 1
50sage: f = cos(x)^6 + sin(x)^6 + 3 * sin(x)^2 * cos(x)^2
51sage: f.simplify_trig()
521
53sage: f = cos(x)^6; f.reduce_trig()
5415/32*cos(2*x) + 3/16*cos(4*x) + 1/32*cos(6*x) + 5/16
55sage: f = sin(5 * x); f.expand_trig()
56sin(x)^5 - 10*sin(x)^3*cos(x)^2 + 5*sin(x)*cos(x)^4
57sage: n = var('n'); f = factorial(n+1)/factorial(n)
58sage: f.simplify_factorial()
59n + 1
60
61"""
62
63
64r"""
65
66sage: reset()
67sage: f = sqrt(x^2); f.simplify_radical()
68abs(x)
69sage: var('y')
70y
71sage: f = log(x*y); f.simplify_radical()
72log(x) + log(y)
73sage: assume(x > 0); bool(sqrt(x^2) == x)
74True
75sage: forget(x > 0); bool(sqrt(x^2) == x)
76False
77sage: var('n'); assume(n, 'integer'); sin(n*pi).simplify()
78n
790
80
81"""
82
83
84r"""
85
86sage: z, phi = var('z, phi')
87sage: solve(z**2 -2 / cos(phi) * z + 5 / cos(phi) ** 2 - 4, z)
88[z == -(2*sqrt(cos(phi)^2 - 1) - 1)/cos(phi), z == (2*sqrt(cos(phi)^2 - 1) + 1)/cos(phi)]
89sage: var('y'); solve(y^6==y, y)
90y
91[y == e^(2/5*I*pi), y == e^(4/5*I*pi), y == e^(-4/5*I*pi), y == e^(-2/5*I*pi), y == 1, y == 0]
92sage: solve(x^2-1, x, solution_dict=True)
93[{x: -1}, {x: 1}]
94sage: solve([x+y == 3, 2*x+2*y == 6],x,y)
95[[x == -r1 + 3, y == r1]]
96sage: solve([cos(x)*sin(x) == 1/2, x+y == 0], x, y)
97[[x == 1/4*pi + pi*z30, y == -1/4*pi - pi*z30]]
98sage: solve(x^2+x-1>0,x)
99[[x < -1/2*sqrt(5) - 1/2], [x > 1/2*sqrt(5) - 1/2]]
100
101"""
102
103
104r"""
105
106sage: x, y, z = var('x, y, z')
107sage: solve([x^2 * y * z == 18, x * y^3 * z == 24,\
108             x * y * z^4 == 6], x, y, z)
109[[x == 3, y == 2, z == 1], [x == (1.33721506733 - 2.68548987407*I), y == (-1.70043427146 + 1.05286432575*I), z == (0.932472229404 - 0.361241666187*I)], [x == (1.33721506733 + 2.68548987407*I), y == (-1.70043427146 - 1.05286432575*I), z == (0.932472229404 + 0.361241666187*I)], [x == (-2.55065140719 - 1.57929648863*I), y == (-0.547325980144 + 1.92365128635*I), z == (-0.982973099684 - 0.183749517817*I)], [x == (-2.55065140719 + 1.57929648863*I), y == (-0.547325980144 - 1.92365128635*I), z == (-0.982973099684 + 0.183749517817*I)], [x == (0.27680507839 - 2.98720252889*I), y == (1.47801783444 - 1.34739128729*I), z == (-0.85021713573 - 0.526432162877*I)], [x == (0.27680507839 + 2.98720252889*I), y == (1.47801783444 + 1.34739128729*I), z == (-0.85021713573 + 0.526432162877*I)], [x == (-0.820988970216 + 2.88547692952*I), y == (-1.20526927276 - 1.59603445456*I), z == (0.0922683594633 - 0.995734176295*I)], [x == (-0.820988970216 - 2.88547692952*I), y == (-1.20526927276 + 1.59603445456*I), z == (0.0922683594633 + 0.995734176295*I)], [x == (-1.80790390914 - 2.39405168184*I), y == (0.891476711553 - 1.79032658271*I), z == (0.739008917221 - 0.673695643647*I)], [x == (-1.80790390914 + 2.39405168184*I), y == (0.891476711553 + 1.79032658271*I), z == (0.739008917221 + 0.673695643647*I)], [x == (2.21702675166 + 2.02108693094*I), y == (1.86494445881 + 0.722483332374*I), z == (-0.273662990072 - 0.961825643173*I)], [x == (2.21702675166 - 2.02108693094*I), y == (1.86494445881 - 0.722483332374*I), z == (-0.273662990072 + 0.961825643173*I)], [x == (2.79741668821 - 1.08372499856*I), y == (-1.96594619937 + 0.367499035633*I), z == (-0.602634636379 - 0.79801722728*I)], [x == (2.79741668821 + 1.08372499856*I), y == (-1.96594619937 - 0.367499035633*I), z == (-0.602634636379 + 0.79801722728*I)], [x == (-2.94891929905 + 0.55124855345*I), y == (0.184536718927 + 1.99146835259*I), z == (0.445738355777 - 0.895163291355*I)], [x == (-2.94891929905 - 0.55124855345*I), y == (0.184536718927 - 1.99146835259*I), z == (0.445738355777 + 0.895163291355*I)]]
110
111"""
112
113
114r"""
115
116sage: expr = sin(x) + sin(2 * x) + sin(3 * x)
117sage: solve(expr, x)
118[sin(3*x) == -sin(2*x) - sin(x)]
119sage: find_root(expr, 0.1, pi)
1202.0943951023931957
121sage: f = expr.simplify_trig(); f
1222*(2*cos(x)^2 + cos(x))*sin(x)
123sage: solve(f, x)
124[x == 0, x == 2/3*pi, x == 1/2*pi]
125sage: (x^3+2*x+1).roots(x)
126[(-1/2*(I*sqrt(3) + 1)*(1/18*sqrt(3)*sqrt(59) - 1/2)^(1/3) - 1/3*(I*sqrt(3) - 1)/(1/18*sqrt(3)*sqrt(59) - 1/2)^(1/3), 1), (-1/2*(-I*sqrt(3) + 1)*(1/18*sqrt(3)*sqrt(59) - 1/2)^(1/3) - 1/3*(-I*sqrt(3) - 1)/(1/18*sqrt(3)*sqrt(59) - 1/2)^(1/3), 1), ((1/18*sqrt(3)*sqrt(59) - 1/2)^(1/3) - 2/3/(1/18*sqrt(3)*sqrt(59) - 1/2)^(1/3), 1)]
127sage: (x^3+2*x+1).roots(x, ring=RR)
128 [(-0.453397651516404, 1)]
129sage: (x^3+2*x+1).roots(x, ring=CC)
130 [(-0.453397651516404, 1), (0.226698825758202 - 1.46771150871022*I, 1), (0.226698825758202 + 1.46771150871022*I, 1)]
131
132"""
133
134
135r"""
136
137sage: var('k n'); sum(k, k, 1, n).factor()
138 (k, n)
139 1/2*(n + 1)*n
140sage: var('n y'); sum(binomial(n,k) * x^k * y^(n-k), k, 0, n)
141 (n, y)
142 (x + y)^n
143sage: var('k n'); sum(binomial(n,k), k, 0, n),\
144                  sum(k * binomial(n, k), k, 0, n),\
145                  sum((-1)^k*binomial(n,k), k, 0, n)
146(k, n)
147(2^n, n*2^(n - 1), 0)
148sage: var('a q'); sum(a*q^k, k, 0, n)
149(a, q)
150(a*q^(n + 1) - a)/(q - 1)
151sage: assume(abs(q) < 1)
152sage: sum(a*q^k, k, 0, infinity)
153 -a/(q - 1)
154
155"""
156
157
158r"""
159
160sage: limit((x**(1/3) -2) / ((x + 19)**(1/3) - 3), x = 8)
1619/4
162sage: f(x) = (cos(pi/4-x)-tan(x))/(1-sin(pi/4 + x))
163sage: limit(f(x), x = pi/4, dir='minus')
164+Infinity
165sage: limit(f(x), x = pi/4, dir='plus')
166-Infinity
167sage: u(n) = n^100 / 100.^n
168sage: u(1);u(2);u(3);u(4);u(5);u(6);u(7);u(8);u(9);u(10)
1690.0100000000000000
1701.26765060022823e26
1715.15377520732011e41
1721.60693804425899e52
1737.88860905221012e59
1746.53318623500071e65
1753.23447650962476e70
1762.03703597633449e74
1772.65613988875875e77
1781.00000000000000e80
179sage: plot(u(x), x, 1, 40)
180sage: v(x) = diff(u(x), x); sol = solve(v(x) == 0, x); sol
181[x == 0, x == (268850/12381)]
182sage: sol[1].rhs().n().floor()
18321
184sage: limit(u(n), n=infinity)
1850
186sage: n0 = find_root(u(n) - 1e-8 == 0, 22, 1000); n0
187105.07496210187252
188
189"""
190
191
192r"""
193
194sage: taylor((1+arctan(x))**(1/x), x, 0, 3)
1951/16*x^3*e + 1/8*x^2*e - 1/2*x*e + e
196sage: (ln(2* sin(x))).series(x==pi/6, 3)
197(sqrt(3))*(-1/6*pi + x) + (-2)*(-1/6*pi + x)^2 + Order(-1/216*(pi - 6*x)^3)
198sage: (ln(2* sin(x))).series(x==pi/6, 3).truncate()
199-1/18*(pi - 6*x)^2 - 1/6*(pi - 6*x)*sqrt(3)
200sage: taylor((x**3+x)**(1/3) - (x**3-x)**(1/3), x, infinity, 2)
2012/3/x
202sage: tan(4*arctan(1/5)).simplify_trig()
203120/119
204sage: tan(pi/4+arctan(1/239)).simplify_trig()
205120/119
206sage: f = arctan(x).series(x, 10); f
2071*x + (-1/3)*x^3 + 1/5*x^5 + (-1/7)*x^7 + 1/9*x^9 + Order(x^10)
208sage: (16*f.subs(x==1/5) - 4*f.subs(x==1/239)).n(); pi.n()
2093.14159268240440
2103.14159265358979
211
212"""
213
214
215r"""
216
217sage: var('k'); sum(1/k^2, k, 1, infinity),\
218                sum(1/k^4, k, 1, infinity),\
219                sum(1/k^5, k, 1, infinity)
220k
221(1/6*pi^2, 1/90*pi^4, zeta(5))
222sage: s = 2*sqrt(2)/9801*(sum((factorial(4*k))*(1103+26390*k)\
223     /((factorial(k)) ^ 4 * 396 ^ (4 * k)) for k in (0..11)))
224sage: (1/s).n(digits=100), pi.n(digits=100), (pi-1/s).n() # abs tol 1e-95
225(3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342121432, 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068, -4.44089209850063e-16)
226sage: var('n'); u = sin(pi*(sqrt(4*n^2+1)-2*n))
227n
228sage: taylor(u, n, infinity, 3)
2291/4*pi/n - 1/384*(6*pi + pi^3)/n^3
230sage: diff(sin(x^2), x)
2312*x*cos(x^2)
232
233"""
234
235
236r"""
237
238sage: function('f',x); function('g',x); diff(f(g(x)), x)
239f(x)
240g(x)
241D[0](f)(g(x))*D[0](g)(x)
242sage: diff(ln(f(x)), x)
243D[0](f)(x)/f(x)
244sage: f(x,y) = x*y + sin(x^2) + e^(-x); derivative(f, x)
245(x, y) |--> 2*x*cos(x^2) + y - e^(-x)
246sage: derivative(f, y)
247(x, y) |--> x
248sage: var('x y'); f = ln(x**2+y**2) / 2
249(x, y)
250sage: delta = diff(f,x,2)+diff(f,y,2)
251sage: delta.simplify_full()
2520
253sage: sin(x).integral(x, 0, pi/2)
2541
255sage: integrate(1/(1+x^2), x)
256arctan(x)
257sage: integrate(1/(1+x^2), x, -infinity, infinity)
258pi
259sage: integrate(exp(-x**2), x, 0, infinity)
2601/2*sqrt(pi)
261sage: u = var('u'); f = x * cos(u) / (u^2 + x^2)
262sage: assume(x>0); f.integrate(u, 0, infinity)
2631/2*pi*e^(-x)
264sage: forget(); assume(x<0); f.integrate(u, 0, infinity)
265-1/2*pi*e^x
266sage: integral_numerical(sin(x)/x, 0, 1)   # abs tol 1e-12
267(0.94608307036718287, 1.0503632079297086e-14)
268sage: g = integrate(exp(-x**2), x, 0, infinity)
269sage: g, g.n()                             # abs tol 1e-12
270(1/2*sqrt(pi), 0.886226925452758)
271sage: approx = integral_numerical(exp(-x**2), 0, infinity)
272sage: approx                               # abs tol 1e-12
273(0.88622692545275705, 1.7147744320162414e-08)
274sage: approx[0]-g.n()                      # abs tol 1e-12
275-8.88178419700125e-16
276
277"""
278
279
280r"""
281
282sage: A = matrix(QQ,[[1,2],[3,4]]); A
283[1 2]
284[3 4]
285sage: A = matrix(QQ, [[2, 4, 3],[-4,-6,-3],[3,3,1]])
286sage: A.characteristic_polynomial()
287x^3 + 3*x^2 - 4
288sage: A.eigenvalues()
289[1, -2, -2]
290sage: A.minimal_polynomial().factor()
291(x - 1) * (x + 2)^2
292sage: A.eigenvectors_right()
293[(1, [
294(1, -1, 1)
295], 1), (-2, [
296(1, -1, 0)
297], 2)]
298sage: A.jordan_form(transformation=True)
299(
300[ 1| 0  0]
301[--+-----]  [ 1  1  1]
302[ 0|-2  1]  [-1 -1  0]
303[ 0| 0 -2], [ 1  0 -1]
304)
305
306"""
307
308
309r"""
310
311sage: A = matrix(QQ, [[1,-1/2],[-1/2,-1]])
312sage: A.minimal_polynomial()
313x^2 - 5/4
314sage: R = QQ[sqrt(5)]
315sage: A = A.change_ring(R)
316sage: A.jordan_form(transformation=True, subdivide=False)
317(
318[ 1/2*sqrt5          0]  [         1          1]
319[         0 -1/2*sqrt5], [-sqrt5 + 2  sqrt5 + 2]
320)
321sage: K.<sqrt2> = NumberField(x^2 - 2)
322sage: L.<sqrt3> = K.extension(x^2 - 3)
323sage: A = matrix(L, [[2, sqrt2*sqrt3, sqrt2], [sqrt2*sqrt3, 3, sqrt3], [sqrt2, sqrt3, 1]])
324sage: A.jordan_form(transformation=True)
325(
326[6|0|0]
327[-+-+-]
328[0|0|0]  [              1               1               0]
329[-+-+-]  [1/2*sqrt2*sqrt3               0               1]
330[0|0|0], [      1/2*sqrt2          -sqrt2          -sqrt3]
331)
332
333"""
334