︠4f58aad8-aa03-4417-90a5-ffaec91c96das︠ ### Manipulating Symbolic Expressions ︡e4e5a99d-0674-401a-ae9b-529bef593415︡{"done":true}︡ ︠407baad3-1d24-44dd-835d-691c4a2d02c5s︠ ### define a symbolic expression f(x) = x+1 type(f) f(x=1) f(1) ︡fccb6aad-488f-48a6-8d72-4010ef9acdb3︡{"stdout":"\n"}︡{"stdout":"2\n"}︡{"stdout":"2\n"}︡{"done":true}︡ ︠9e597852-aa2e-4bc5-9a45-3d9a67a8c0e7︠ ### Better to tell the name of the variable f = x+1 f(1) ︡9232e7af-af4f-4384-9294-0bc967e848e7︡{"stdout":"2\n"}︡{"stderr":"/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py:1188: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...)\nSee http://trac.sagemath.org/5930 for details.\n flags=compile_flags) in namespace, locals\n"}︡{"done":true}︡ ︠eea2a0e4-d7a3-4947-baf5-aba3b494d3b4s︠ ### x is a pre-defined variable while y is not f = x+y ︡e1e2a099-f0b4-4639-acf1-63bf0bc66376︡{"stderr":"Error in lines 1-1\nTraceback (most recent call last):\n File \"/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py\", line 1188, in execute\n flags=compile_flags) in namespace, locals\n File \"\", line 1, in \nNameError: name 'y' is not defined\n"}︡{"done":true}︡ ︠81557aaa-0fbd-45d3-b4a3-b8ce7a126c7ds︠ ### set up some variables y = var('y') ### create a variable named 'y', and store it into a Python variable y f = x+y ︡8203d2be-f39c-4b83-943c-9c6b02c0ee10︡{"done":true}︡ ︠8053c892-4ec1-458a-8ee1-6a7f2f3f137f︠ f(1) ### default will substitute x ︡dd71ec19-3b50-4c9e-84d0-756f7bf4c9ba︡{"stdout":"y + 1\n"}︡{"done":true}︡ ︠1c1bc840-2a8c-475d-868a-4381918571ees︠ f(y=1) ### but you can specifically tell which one to substitute ︡026b661a-f13c-4915-b69e-f281586439fb︡{"stdout":"x + 1\n"}︡{"done":true}︡ ︠368180d7-68da-4429-a7c3-8ddac8825416s︠ ### differentiate f = x^2 + y^3 f.differentiate(x) f.differentiate(y) ︡8e0da1c0-0bef-449a-b504-67b59dc25e09︡{"stdout":"2*x\n"}︡{"stdout":"3*y^2\n"}︡{"done":true}︡ ︠a4748a0d-3f6d-4574-addc-304d994e1675s︠ ### rational functions f = (x^10-1)/(x-1) f(1) ### Sage will complain when the divisor is 0 ︡4ff816cc-488b-4601-bd1d-cdd18b136c89︡{"stderr":"Error in lines 2-2\nTraceback (most recent call last):\n File \"/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py\", line 1188, in execute\n flags=compile_flags) in namespace, locals\n File \"\", line 1, in \n File \"sage/symbolic/expression.pyx\", line 5467, in sage.symbolic.expression.Expression.__call__ (build/cythonized/sage/symbolic/expression.cpp:34053)\n return self._parent._call_element_(self, *args, **kwds)\n File \"sage/symbolic/ring.pyx\", line 980, in sage.symbolic.ring.SymbolicRing._call_element_ (build/cythonized/sage/symbolic/ring.cpp:11750)\n return _the_element.subs(d, **kwds)\n File \"sage/symbolic/expression.pyx\", line 5319, in sage.symbolic.expression.Expression.substitute (build/cythonized/sage/symbolic/expression.cpp:32983)\n res = self._gobj.subs_map(smap, 0)\nValueError: power::eval(): division by zero\n"}︡{"done":true}︡ ︠22c034de-2ac7-40c2-811b-c62989e90816︠ ### Numerical solution: shed light on limit for k in [1.1,1.01,1.001,1.0001]: print k, f(x=k); ︡7a2815c0-b8bd-4a6b-84ca-1037a75354ad︡{"stdout":"1.10000000000000 15.9374246010000\n1.01000000000000 10.4622125411204\n1.00100000000000 10.0451202102523\n1.00010000000000 10.0045012002103\n"}︡{"done":true}︡ ︠169fcbc4-3d96-4300-b448-ad402552e5b4s︠ ### Alternative solution: L'Hopital's rule ### thinking f = p / q p = f.numerator(normalize=False) q = f.denominator(normalize=False) print "p =", p print "p(1) =", p(x=1) print "q =", q print "q(1) =", q(x=1) ︡9e6b3178-2473-4bf7-ab9b-327969a3fc58︡{"stdout":"p = (x^10 - 1)\n"}︡{"stdout":"p(1) = 0\n"}︡{"stdout":"q = x - 1\n"}︡{"stdout":"q(1) = 0\n"}︡{"done":true}︡ ︠ed207632-ed47-44fb-8575-440db5b8a466︠ ### What does normalize do? p = f.numerator(normalize=True) q = f.denominator(normalize=True) print "p =", p print "p(1) =", p(x=1) print "q =", q print "q(1) =", q(x=1) ︡1cf4ea1f-af6d-4b93-9b00-b46e1c6873f7︡{"stdout":"p = x^9 + x^8 + x^7 + x^6 + x^5 + x^4 + x^3 + x^2 + x + 1\n"}︡{"stdout":"p(1) = 10\n"}︡{"stdout":"q = 1\n"}︡{"stdout":"q(1) = 1\n"}︡{"done":true}︡ ︠870034df-40d3-4992-bf27-44ff8298231as︠ ### L'Hopital's rule ### thinking f = p / q p = f.numerator(normalize=False) q = f.denominator(normalize=False) print "p =", p print "p(1) =", p(x=1) print "q =", q print "q(1) =", q(x=1) print "---" pprime = p.differentiate(x) qprime = q.differentiate(x) print "p' =", pprime print "q' =", qprime print "p'(1) =", pprime(x=1) print "q'(1) =", qprime(x=1) print "---" print "so the limit is {}".format(pprime(x=1)/qprime(x=1)) ︡2aa2da92-7dba-450d-8636-347ff45e2ddb︡{"stdout":"p = (x^10 - 1)\n"}︡{"stdout":"p(1) = 0\n"}︡{"stdout":"q = x - 1\n"}︡{"stdout":"q(1) = 0\n"}︡{"stdout":"---\n"}︡{"stdout":"p' = 10*x^9\n"}︡{"stdout":"q' = 1\n"}︡{"stdout":"p'(1) = 10\n"}︡{"stdout":"q'(1) = 1\n"}︡{"stdout":"---\n"}︡{"stdout":"so the limit is 10\n"}︡{"done":true}︡ ︠1534de01-8b7d-4730-b62d-82ed3fb994f6s︠ ### Let's do the tedious work ### goal: let f = (xD)^2 (x^(N+1) - 1)/(x - 1) ### find lim f when x approaches 1. N = var('N') x = var('x') b = (x^(N+1) - 1) / (x - 1) b(x=1) ︡12cda28a-6777-4135-b409-097a6f0c943a︡{"stderr":"Error in lines 4-4\nTraceback (most recent call last):\n File \"/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py\", line 1188, in execute\n flags=compile_flags) in namespace, locals\n File \"\", line 1, in \n File \"sage/symbolic/expression.pyx\", line 5467, in sage.symbolic.expression.Expression.__call__ (build/cythonized/sage/symbolic/expression.cpp:34053)\n return self._parent._call_element_(self, *args, **kwds)\n File \"sage/symbolic/ring.pyx\", line 980, in sage.symbolic.ring.SymbolicRing._call_element_ (build/cythonized/sage/symbolic/ring.cpp:11750)\n return _the_element.subs(d, **kwds)\n File \"sage/symbolic/expression.pyx\", line 5319, in sage.symbolic.expression.Expression.substitute (build/cythonized/sage/symbolic/expression.cpp:32983)\n res = self._gobj.subs_map(smap, 0)\nValueError: power::eval(): division by zero\n"}︡{"done":true}︡ ︠49a49323-2614-45da-a902-3c93fedeb61cs︠ ### not good for LH rule b.differentiate(x).show() ︡f197189c-5436-4a4c-891d-9f59c3fa2d44︡{"html":"
$\\displaystyle \\frac{{\\left(N + 1\\right)} x^{N}}{x - 1} - \\frac{x^{N + 1} - 1}{{\\left(x - 1\\right)}^{2}}$
"}︡{"done":true}︡ ︠9d476447-04c4-45d6-8a3b-99a7b230d1cfs︠ ### normalize combines the denominators: good for LH rule b.differentiate(x).normalize().show() ︡c9d2d659-97f8-4f47-9383-e67370a2d41f︡{"html":"
$\\displaystyle \\frac{N x x^{N} - N x^{N} + x x^{N} - x^{N + 1} - x^{N} + 1}{{\\left(x - 1\\right)}^{2}}$
"}︡{"done":true}︡ ︠0106e4cd-24fd-4f99-a4bd-87502d10ac55s︠ ### the numerators are not simplified ### define the function below to simplify the denominator def clean_numerator(p): return p.numerator().expand().simplify() / p.denominator() ︡e173814e-e349-4365-952e-df17e9fb8673︡{"done":true}︡ ︠fcb700d5-883e-4bb0-b348-f383a85bf87fs︠ b1 = x * b.differentiate(x).normalize() b1.show() ︡2c800432-89b9-4983-b78f-ee2ccad9d7d6︡{"html":"
$\\displaystyle \\frac{{\\left(N x x^{N} - N x^{N} + x x^{N} - x^{N + 1} - x^{N} + 1\\right)} x}{{\\left(x - 1\\right)}^{2}}$
"}︡{"done":true}︡ ︠007f36e7-8bf8-4441-9e2e-fa96ff5084e4s︠ b1 = clean_numerator(b1) b1.show() ︡02b6d780-c4c3-4961-b9f2-d50588096e49︡{"html":"
$\\displaystyle \\frac{N x^{N + 2} - N x^{N + 1} + x - x^{N + 1}}{{\\left(x - 1\\right)}^{2}}$
"}︡{"done":true}︡ ︠e0837ac3-dbdc-470d-94b0-c39c497e8801s︠ b2 = x * b1.differentiate(x).normalize() b2.show() ︡8e058d2a-25bf-4cc2-8d10-de788fa9939c︡{"html":"
$\\displaystyle \\frac{{\\left(N^{2} x x^{N + 1} - N^{2} x x^{N} - N^{2} x^{N + 1} + 2 \\, N x x^{N + 1} + N^{2} x^{N} - 2 \\, N x x^{N} - 2 \\, N x^{N + 2} + 2 \\, N x^{N} - x x^{N} + 2 \\, x^{N + 1} + x^{N} - x - 1\\right)} x}{{\\left(x - 1\\right)}^{3}}$
"}︡{"done":true}︡ ︠5859bfbd-2325-449e-8317-7d27868f01f6s︠ b2 = clean_numerator(b2) b2.show() ︡222109c0-b273-42df-b401-6a0b28c60f5d︡{"html":"
$\\displaystyle \\frac{N^{2} x^{N + 3} - 2 \\, N^{2} x^{N + 2} + N^{2} x^{N + 1} - 2 \\, N x^{N + 2} + 2 \\, N x^{N + 1} - x^{2} + x^{N + 2} + x^{N + 1} - x}{{\\left(x - 1\\right)}^{3}}$
"}︡{"done":true}︡ ︠55fea383-17b4-4dad-8541-349e1c98ea2cs︠ num = b2.numerator() den = b2.denominator() print "num =", num print "num(1) =", num(x=1) print "den =", den print "den(1) =", den(x=1) ︡aed41083-a9b5-4aef-a1ef-7a8345d153e9︡{"stdout":"num = N^2*x^(N + 3) - 2*N^2*x^(N + 2) + N^2*x^(N + 1) - 2*N*x^(N + 2) + 2*N*x^(N + 1) - x^2 + x^(N + 2) + x^(N + 1) - x\n"}︡{"stdout":"num(1) = 0\n"}︡{"stdout":"den = (x - 1)^3\n"}︡{"stdout":"den(1) = 0\n"}︡{"done":true}︡ ︠2d95eea2-cacf-415d-a316-73e93cfcdf70s︠ ### Apply LH rule 3 times numppp = num.differentiate(x,3) denppp = den.differentiate(x,3) print "num''' =", numppp print "num'''(1) =", numppp(x=1) print "den''' =", denppp print "den'''(1) =", denppp(x=1) ︡6a5a3539-b6e7-4646-9167-d7469fa56177︡{"stdout":"num''' = -2*(N + 2)*(N + 1)*N^3*x^(N - 1) + (N + 1)*(N - 1)*N^3*x^(N - 2) + (N + 3)*(N + 2)*(N + 1)*N^2*x^N - 2*(N + 2)*(N + 1)*N^2*x^(N - 1) + 2*(N + 1)*(N - 1)*N^2*x^(N - 2) + (N + 2)*(N + 1)*N*x^(N - 1) + (N + 1)*(N - 1)*N*x^(N - 2)\n"}︡{"stdout":"num'''(1) = (N + 3)*(N + 2)*(N + 1)*N^2 - 2*(N + 2)*(N + 1)*N^3 + (N + 1)*(N - 1)*N^3 - 2*(N + 2)*(N + 1)*N^2 + 2*(N + 1)*(N - 1)*N^2 + (N + 2)*(N + 1)*N + (N + 1)*(N - 1)*N\n"}︡{"stdout":"den''' = 6\n"}︡{"stdout":"den'''(1) = 6\n"}︡{"done":true}︡ ︠cf1d14f7-5adb-4347-9d26-c1ceb17bd0b0s︠ last = numppp(x=1).expand().simplify() last.show() ︡458771e6-6e56-49c8-a490-dc3a0fdfee71︡{"html":"
$\\displaystyle 2 \\, N^{3} + 3 \\, N^{2} + N$
"}︡{"done":true}︡ ︠1269c789-083e-435f-a386-46baabfc148a︠ last.factor().show() ### Bingo!!! ︡cf4aad4a-5b93-4f21-9579-763f00628af6︡{"html":"
$\\displaystyle {\\left(2 \\, N + 1\\right)} {\\left(N + 1\\right)} N$
"}︡{"done":true}︡ ︠9f20d4a5-6626-47fd-9308-d8da5f58503di︠ %md ## Project 1 Following the same idea, find a formula for $\sum_{k=1}^N k^3$. Bonus: Can you write a function `power_sum(p)` to compute the formula of $\sum_{k=1}^N k^p$. ︡1491331a-ff0e-4231-ae95-9c18acf6be1f︡{"done":true,"md":"\n## Project 1\n\nFollowing the same idea, find a formula for $\\sum_{k=1}^N k^3$. \nBonus: Can you write a function `power_sum(p)` to compute the formula of $\\sum_{k=1}^N k^p$."} ︠49a8da51-288f-4181-b23f-c90972e0f3f6s︠ def clean_numerator(p): return p.numerator().expand().simplify() / p.denominator() def xD(f): return clean_numerator(x * f.differentiate(x).normalize()) def power_sum(p): ### Your answer here N = var('N') x = var('x') b = (x^(N+1) - 1) / (x - 1) for _ in range(p): b = xD(b) num = b.numerator() den = b.denominator() last_num = num.differentiate(x,p+1).subs(x=1) last_den = den.differentiate(x,p+1).subs(x=1) return last_num.expand().simplify().factor() / last_den ︡16896951-8517-41e3-93e9-86be40679397︡{"done":true}︡ ︠948c4818-f37f-4852-b6e6-2f9b15b0897ds︠ for p in range(5): print power_sum(p) ︡6c3282aa-bd9f-44dc-baad-425f69e5c022︡{"stdout":"N + 1\n1/2*(N + 1)*N\n1/6*(2*N + 1)*(N + 1)*N\n1/4*(N + 1)^2*N^2\n1/30*(3*N^2 + 3*N - 1)*(2*N + 1)*(N + 1)*N\n"}︡{"done":true}︡ ︠65a5278c-874d-4f3b-942e-af546124bdec︠ ### Some answers for references ( p = 0,1,2,3,4) N + 1 1/2*(N + 1)*N 1/6*(2*N + 1)*(N + 1)*N 1/4*(N + 1)^2*N^2 1/30*(3*N^2 + 3*N - 1)*(2*N + 1)*(N + 1)*N