Skip to content

Instantly share code, notes, and snippets.

@amca01

amca01/solera.ipynb

Created December 12, 2023 00:44
Show Gist options
  • Save amca01/362ac302882446bf2edffb7ee5efc3e4 to your computer and use it in GitHub Desktop.
Save amca01/362ac302882446bf2edffb7ee5efc3e4 to your computer and use it in GitHub Desktop.
Download ZIP
Solera jupyter notebook
Display the source blob
Display the rendered blob
Raw
solera.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 76,
"id": "564fc7f9-fccc-4b29-94e0-74e64f17fb42",
"metadata": {},
"outputs": [],
"source": [
"import sympy as sy"
]
},
{
"cell_type": "code",
"execution_count": 97,
"id": "5758370d-f10c-4356-996d-85dd8118b79c",
"metadata": {},
"outputs": [],
"source": [
"a,b,c,d,e,f,g = sy.symbols('a,b,c,d,e,f,g',cls=sy.Function)\n",
"n,p = sy.symbols('n,p')"
]
},
{
"cell_type": "code",
"execution_count": 98,
"id": "b00b0a7a-5089-4be1-b488-7c2bc25866c0",
"metadata": {},
"outputs": [],
"source": [
"def R(x,y):\n",
" return sy.Rational(x,y)"
]
},
{
"cell_type": "markdown",
"id": "1540516b-8174-464f-8b5e-ed51767f32c7",
"metadata": {},
"source": [
"# Numerical experiments"
]
},
{
"cell_type": "code",
"execution_count": 143,
"id": "773a94fc-09ed-4f52-a816-8cbb69f800aa",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}1 & 2.99999999999864 & 5.99999927650091 & 9.99970047294786\\\\1 & 2.99999999999932 & 5.99999951766682 & 9.99977517383613\\\\1 & 2.99999999999966 & 5.99999967844432 & 9.9998312597938\\\\1 & 2.99999999999983 & 5.99999978562943 & 9.99987336445643\\\\1 & 2.99999999999991 & 5.99999985708623 & 9.99990496974968\\\\1 & 2.99999999999996 & 5.99999990472413 & 9.99992869158382\\\\1 & 2.99999999999998 & 5.99999993648274 & 9.9999464948689\\\\1 & 2.99999999999999 & 5.99999995765515 & 9.99995985527236\\\\1 & 2.99999999999999 & 5.9999999717701 & 9.99996988086805\\\\1 & 3.0 & 5.99999998118006 & 9.99997740359357\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[1, 2.99999999999864, 5.99999927650091, 9.99970047294786],\n",
"[1, 2.99999999999932, 5.99999951766682, 9.99977517383613],\n",
"[1, 2.99999999999966, 5.99999967844432, 9.9998312597938],\n",
"[1, 2.99999999999983, 5.99999978562943, 9.99987336445643],\n",
"[1, 2.99999999999991, 5.99999985708623, 9.99990496974968],\n",
"[1, 2.99999999999996, 5.99999990472413, 9.99992869158382],\n",
"[1, 2.99999999999998, 5.99999993648274, 9.9999464948689],\n",
"[1, 2.99999999999999, 5.99999995765515, 9.99995985527236],\n",
"[1, 2.99999999999999, 5.9999999717701, 9.99996988086805],\n",
"[1, 3.0, 5.99999998118006, 9.99997740359357]])"
]
},
"execution_count": 143,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"abcd = sy.Matrix([[1,0,0,0]])\n",
"for i in range(50):\n",
" td = 3/4*abcd[-1,3]+1/4*abcd[-1,2]+1\n",
" tc = 2/3*abcd[-1,2]+1/3*abcd[-1,1]+1\n",
" tb = 1/2*abcd[-1,1]+3/2\n",
" ta = 1\n",
" tr = sy.Matrix([[ta,tb,tc,td]])\n",
" abcd = abcd.col_join(tr)\n",
"abcd[-10:,:]\n",
" "
]
},
{
"cell_type": "markdown",
"id": "5345ceda-3efc-4761-9bcc-73fa3acccee8",
"metadata": {},
"source": [
"## Solving the equations using \"rsolve\""
]
},
{
"cell_type": "code",
"execution_count": 99,
"id": "f2a514f1-4d1e-4239-b3c2-fcb4d361654b",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 3 - 3 \\cdot 2^{- n}$"
],
"text/plain": [
"3 - 3/2**n"
]
},
"execution_count": 99,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bn = sy.rsolve(2*b(n+1)-b(n)-3,b(n),{b(0):0})\n",
"bn"
]
},
{
"cell_type": "code",
"execution_count": 105,
"id": "1482eb11-ed09-4284-8c6c-83010ae3c994",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 6 - 12 \\left(\\frac{3}{2}\\right)^{- n} + 6 \\cdot 2^{- n}$"
],
"text/plain": [
"6 - 12/(3/2)**n + 6/2**n"
]
},
"execution_count": 105,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cn = sy.rsolve(3*c(n+1)-2*c(n)-bn-3,c(n),{c(0):0})\n",
"cn"
]
},
{
"cell_type": "code",
"execution_count": 115,
"id": "232d9ec5-a033-4652-ae49-a6afd55b601d",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 10 - 40 \\left(\\frac{4}{3}\\right)^{- n} + 36 \\left(\\frac{3}{2}\\right)^{- n} - 6 \\cdot 2^{- n}$"
],
"text/plain": [
"10 - 40/(4/3)**n + 36/(3/2)**n - 6/2**n"
]
},
"execution_count": 115,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dn = sy.rsolve(4*d(n+1)-3*d(n)-cn-4,d(n),{d(0):0})\n",
"dn.expand().powsimp(combine='base')#,force=True)"
]
},
{
"cell_type": "markdown",
"id": "bc292841-2367-42de-be5d-e82608f9460b",
"metadata": {},
"source": [
"# Creating a single difference equation for each row\n",
"\n",
"From\n",
"\n",
"(1) $b_{n+1}=\\frac{1}{2}b_n+3/2$\n",
"\n",
"(2) $c_{n+1}=\\frac{2}{3}c_n+\\frac{1}{3}b_n+1$\n",
"\n",
"(3) $d_{n+1}=\\frac{3}{4}d_n+\\frac{1}{4}c_n+1$\n",
"\n",
"From (2):\n",
"\n",
"(2') $b_n=3c_{n+1}-2c_n-3$\n",
"\n",
"Substitute back into (1):\n",
"\n",
"$3c_{n+2}-2c_{n+1}-3=\\frac{1}{2}(3c_{n+1}-2c_n-3)+3/2$\n",
"\n",
"or without fractions:\n",
"\n",
"$6c_{n+2}-4c_{n+1}-6=3c_{n+1}-2c_n-3+3\\quad\\Rightarrow\\quad 6c_{n+2}-7c_{n+1}+2c_n=6.$"
]
},
{
"cell_type": "code",
"execution_count": 116,
"id": "c1c87045-951a-4367-bb0f-ac6a0ba2e6e8",
"metadata": {},
"outputs": [],
"source": [
"b0,b1,c0,c1,c2,d0,d1,d2,d3= sy.symbols('b0,b1,c0,c1,c2,d0,d1,d2,d3')"
]
},
{
"cell_type": "code",
"execution_count": 122,
"id": "c6725219-4125-4337-b8d6-630bdbc4d7ff",
"metadata": {},
"outputs": [],
"source": [
"bc = {b0:3*c1-2*c0-3,b1:3*c2-2*c1-3}"
]
},
{
"cell_type": "code",
"execution_count": 133,
"id": "d32e645b-4c21-4f7b-85e8-1d81ff6a467b",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 2 c_{0} - 7 c_{1} + 6 c_{2} - 6$"
],
"text/plain": [
"2*c0 - 7*c1 + 6*c2 - 6"
]
},
"execution_count": 133,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cn1 = (2*b1-b0-3).subs(bc)\n",
"cn1"
]
},
{
"cell_type": "code",
"execution_count": 130,
"id": "c9217522-af7d-48ea-b832-41fffbc82a8d",
"metadata": {},
"outputs": [],
"source": [
"cd = {c0:4*d1-3*d0-4,c1:4*d2-3*d1-4,c2:4*d3-3*d2-4}"
]
},
{
"cell_type": "code",
"execution_count": 131,
"id": "208502d6-8074-411c-80ed-e3f24d1e99a2",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle - 6 d_{0} + 29 d_{1} - 46 d_{2} + 24 d_{3} - 10$"
],
"text/plain": [
"-6*d0 + 29*d1 - 46*d2 + 24*d3 - 10"
]
},
"execution_count": 131,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cn1.subs(cd)"
]
},
{
"cell_type": "code",
"execution_count": 182,
"id": "9eebfdf2-07a6-48ac-ab57-5e16bbdc5a05",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}0 & 0 & 0 & 0 & 0 & 0 & 0\\\\1 & 0 & 0 & 0 & 0 & 0 & 0\\\\2 & -1 & 0 & 0 & 0 & 0 & 0\\\\6 & -7 & 2 & 0 & 0 & 0 & 0\\\\24 & -46 & 29 & -6 & 0 & 0 & 0\\\\120 & -326 & 329 & -146 & 24 & 0 & 0\\\\720 & -2556 & 3604 & -2521 & 874 & -120 & 0\\\\5040 & -22212 & 40564 & -39271 & 21244 & -6084 & 720\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ 0, 0, 0, 0, 0, 0, 0],\n",
"[ 1, 0, 0, 0, 0, 0, 0],\n",
"[ 2, -1, 0, 0, 0, 0, 0],\n",
"[ 6, -7, 2, 0, 0, 0, 0],\n",
"[ 24, -46, 29, -6, 0, 0, 0],\n",
"[ 120, -326, 329, -146, 24, 0, 0],\n",
"[ 720, -2556, 3604, -2521, 874, -120, 0],\n",
"[5040, -22212, 40564, -39271, 21244, -6084, 720]])"
]
},
"execution_count": 182,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k = 6\n",
"T = sy.zeros(k+2,k+1)\n",
"T[1,0]=1\n",
"for i in range(2,k+2):\n",
" T[i,0] = i*T[i-1,0]\n",
" for j in range(1,k+1):\n",
" T[i,j] = i*T[i-1,j]-(i-1)*T[i-1,j-1]\n",
"T"
]
},
{
"cell_type": "code",
"execution_count": 183,
"id": "aa0369cd-3723-4047-a2e6-66000a56f158",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}0\\\\1\\\\1\\\\1\\\\1\\\\1\\\\1\\\\1\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[0],\n",
"[1],\n",
"[1],\n",
"[1],\n",
"[1],\n",
"[1],\n",
"[1],\n",
"[1]])"
]
},
"execution_count": 183,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"T*sy.ones(k+1,1)"
]
},
{
"cell_type": "code",
"execution_count": 159,
"id": "2b7f10c7-0b37-49b5-804e-00e65af8ed2e",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 2 x - 1$"
],
"text/plain": [
"2*x - 1"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle 2 x - 1$"
],
"text/plain": [
"2*x - 1"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[1/2]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle 6 x^{2} - 7 x + 2$"
],
"text/plain": [
"6*x**2 - 7*x + 2"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\left(2 x - 1\\right) \\left(3 x - 2\\right)$"
],
"text/plain": [
"(2*x - 1)*(3*x - 2)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[1/2, 2/3]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle 24 x^{3} - 46 x^{2} + 29 x - 6$"
],
"text/plain": [
"24*x**3 - 46*x**2 + 29*x - 6"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\left(2 x - 1\\right) \\left(3 x - 2\\right) \\left(4 x - 3\\right)$"
],
"text/plain": [
"(2*x - 1)*(3*x - 2)*(4*x - 3)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[1/2, 2/3, 3/4]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle 120 x^{4} - 326 x^{3} + 329 x^{2} - 146 x + 24$"
],
"text/plain": [
"120*x**4 - 326*x**3 + 329*x**2 - 146*x + 24"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\left(2 x - 1\\right) \\left(3 x - 2\\right) \\left(4 x - 3\\right) \\left(5 x - 4\\right)$"
],
"text/plain": [
"(2*x - 1)*(3*x - 2)*(4*x - 3)*(5*x - 4)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[1/2, 2/3, 3/4, 4/5]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle 720 x^{5} - 2556 x^{4} + 3604 x^{3} - 2521 x^{2} + 874 x - 120$"
],
"text/plain": [
"720*x**5 - 2556*x**4 + 3604*x**3 - 2521*x**2 + 874*x - 120"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\left(2 x - 1\\right) \\left(3 x - 2\\right) \\left(4 x - 3\\right) \\left(5 x - 4\\right) \\left(6 x - 5\\right)$"
],
"text/plain": [
"(2*x - 1)*(3*x - 2)*(4*x - 3)*(5*x - 4)*(6*x - 5)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[1/2, 2/3, 3/4, 4/5, 5/6]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle 5040 x^{6} - 22212 x^{5} + 40564 x^{4} - 39271 x^{3} + 21244 x^{2} - 6084 x + 720$"
],
"text/plain": [
"5040*x**6 - 22212*x**5 + 40564*x**4 - 39271*x**3 + 21244*x**2 - 6084*x + 720"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\left(2 x - 1\\right) \\left(3 x - 2\\right) \\left(4 x - 3\\right) \\left(5 x - 4\\right) \\left(6 x - 5\\right) \\left(7 x - 6\\right)$"
],
"text/plain": [
"(2*x - 1)*(3*x - 2)*(4*x - 3)*(5*x - 4)*(6*x - 5)*(7*x - 6)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[1/2, 2/3, 3/4, 4/5, 5/6, 6/7]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x = sy.symbols('x')\n",
"for i in range(2,k+2):\n",
" px = sum(T[i,j]*x**(i-j-1) for j in range(i))\n",
" display(px,px.factor(),sy.solve(px,x))"
]
},
{
"cell_type": "markdown",
"id": "9bb562a4-99ad-4070-808e-a8d5af38ae5e",
"metadata": {},
"source": [
"# Using matrices"
]
},
{
"cell_type": "code",
"execution_count": 79,
"id": "6a720df1-4d54-4d43-9ca6-a801ed5a6a9e",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}\\frac{1}{2} & 0 & 0 & 0 & 0\\\\\\frac{1}{3} & \\frac{2}{3} & 0 & 0 & 0\\\\0 & \\frac{1}{4} & \\frac{3}{4} & 0 & 0\\\\0 & 0 & \\frac{1}{5} & \\frac{4}{5} & 0\\\\0 & 0 & 0 & \\frac{1}{6} & \\frac{5}{6}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[1/2, 0, 0, 0, 0],\n",
"[1/3, 2/3, 0, 0, 0],\n",
"[ 0, 1/4, 3/4, 0, 0],\n",
"[ 0, 0, 1/5, 4/5, 0],\n",
"[ 0, 0, 0, 1/6, 5/6]])"
]
},
"execution_count": 79,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = sy.Matrix([[R(1,2),0,0,0,0],[R(1,3),R(2,3),0,0,0],[0,R(1,4),R(3,4),0,0],[0,0,R(1,5),R(4,5),0],[0,0,0,R(1,6),R(5,6)]])\n",
"A"
]
},
{
"cell_type": "code",
"execution_count": 80,
"id": "ea915402-7191-414a-85b2-e589184c46b1",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(Matrix([\n",
" [ 3, 0, 0, 0, 0],\n",
" [-6, -2, 0, 0, 0],\n",
" [ 6, 6, 1, 0, 0],\n",
" [-4, -9, -4, -1, 0],\n",
" [ 2, 9, 8, 5, 1]]),\n",
" Matrix([\n",
" [1/2, 0, 0, 0, 0],\n",
" [ 0, 2/3, 0, 0, 0],\n",
" [ 0, 0, 3/4, 0, 0],\n",
" [ 0, 0, 0, 4/5, 0],\n",
" [ 0, 0, 0, 0, 5/6]]))"
]
},
"execution_count": 80,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P,D = A.diagonalize()\n",
"P,D"
]
},
{
"cell_type": "code",
"execution_count": 81,
"id": "8978e0a3-4c7b-4b02-8d02-18c329974058",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}\\frac{1}{2} & 0 & 0 & 0 & 0\\\\\\frac{1}{3} & \\frac{2}{3} & 0 & 0 & 0\\\\0 & \\frac{1}{4} & \\frac{3}{4} & 0 & 0\\\\0 & 0 & \\frac{1}{5} & \\frac{4}{5} & 0\\\\0 & 0 & 0 & \\frac{1}{6} & \\frac{5}{6}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[1/2, 0, 0, 0, 0],\n",
"[1/3, 2/3, 0, 0, 0],\n",
"[ 0, 1/4, 3/4, 0, 0],\n",
"[ 0, 0, 1/5, 4/5, 0],\n",
"[ 0, 0, 0, 1/6, 5/6]])"
]
},
"execution_count": 81,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P*D*P.inv()"
]
},
{
"cell_type": "code",
"execution_count": 82,
"id": "9ba84c45-9e48-4833-9a38-de5ba13a3cea",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}-2 & 0 & 0 & 0 & 0\\\\-2 & -3 & 0 & 0 & 0\\\\-2 & -3 & -4 & 0 & 0\\\\-2 & -3 & -4 & -5 & 0\\\\-2 & -3 & -4 & -5 & -6\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[-2, 0, 0, 0, 0],\n",
"[-2, -3, 0, 0, 0],\n",
"[-2, -3, -4, 0, 0],\n",
"[-2, -3, -4, -5, 0],\n",
"[-2, -3, -4, -5, -6]])"
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"I5 = sy.eye(5)\n",
"(A - I5).inv()"
]
},
{
"cell_type": "code",
"execution_count": 86,
"id": "4422e3c2-195a-457d-9f45-a6f804c31965",
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}3 - 3 \\cdot 2^{- n}\\\\6 - 12 \\left(\\frac{3}{2}\\right)^{- n} + 6 \\cdot 2^{- n}\\\\10 - 40 \\left(\\frac{4}{3}\\right)^{- n} + 36 \\left(\\frac{3}{2}\\right)^{- n} - 6 \\cdot 2^{- n}\\\\15 - 125 \\left(\\frac{5}{4}\\right)^{- n} + 160 \\left(\\frac{4}{3}\\right)^{- n} - 54 \\left(\\frac{3}{2}\\right)^{- n} + 4 \\cdot 2^{- n}\\\\21 - 378 \\left(\\frac{6}{5}\\right)^{- n} + 625 \\left(\\frac{5}{4}\\right)^{- n} - 320 \\left(\\frac{4}{3}\\right)^{- n} + 54 \\left(\\frac{3}{2}\\right)^{- n} - 2 \\cdot 2^{- n}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ 3 - 3/2**n],\n",
"[ 6 - 12/(3/2)**n + 6/2**n],\n",
"[ 10 - 40/(4/3)**n + 36/(3/2)**n - 6/2**n],\n",
"[ 15 - 125/(5/4)**n + 160/(4/3)**n - 54/(3/2)**n + 4/2**n],\n",
"[21 - 378/(6/5)**n + 625/(5/4)**n - 320/(4/3)**n + 54/(3/2)**n - 2/2**n]])"
]
},
"execution_count": 86,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"An = A**n\n",
"b = sy.Matrix([R(3,2),1,1,1,1])\n",
"(An - I5)*(A - I5).inv()*b"
]
},
{
"cell_type": "code",
"execution_count": 87,
"id": "d22bcbb0-2151-441e-ac99-46dbfbbf0510",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}1 - p & 0 & 0 & 0 & 0 & 0\\\\\\frac{p}{2} & 1 - \\frac{p}{2} & 0 & 0 & 0 & 0\\\\0 & \\frac{p}{3} & 1 - \\frac{p}{3} & 0 & 0 & 0\\\\0 & 0 & \\frac{p}{4} & 1 - \\frac{p}{4} & 0 & 0\\\\0 & 0 & 0 & \\frac{p}{5} & 1 - \\frac{p}{5} & 0\\\\0 & 0 & 0 & 0 & \\frac{p}{6} & 1 - \\frac{p}{6}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[1 - p, 0, 0, 0, 0, 0],\n",
"[ p/2, 1 - p/2, 0, 0, 0, 0],\n",
"[ 0, p/3, 1 - p/3, 0, 0, 0],\n",
"[ 0, 0, p/4, 1 - p/4, 0, 0],\n",
"[ 0, 0, 0, p/5, 1 - p/5, 0],\n",
"[ 0, 0, 0, 0, p/6, 1 - p/6]])"
]
},
"execution_count": 87,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Ap = sy.Matrix([[1-p,0,0,0,0,0],[p/2,1-p/2,0,0,0,0],[0,p/3,1-p/3,0,0,0],[0,0,p/4,1-p/4,0,0],[0,0,0,p/5,1-p/5,0],[0,0,0,0,p/6,1-p/6]])\n",
"Ap"
]
},
{
"cell_type": "code",
"execution_count": 88,
"id": "b8bc0067-8398-4153-bc42-8c9aff6b70f5",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}1 - p & 0 & 0 & 0\\\\\\frac{p}{2} & 1 - \\frac{p}{2} & 0 & 0\\\\0 & \\frac{p}{3} & 1 - \\frac{p}{3} & 0\\\\0 & 0 & \\frac{p}{4} & 1 - \\frac{p}{4}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[1 - p, 0, 0, 0],\n",
"[ p/2, 1 - p/2, 0, 0],\n",
"[ 0, p/3, 1 - p/3, 0],\n",
"[ 0, 0, p/4, 1 - p/4]])"
]
},
"execution_count": 88,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A4 = Ap[:4,:4]\n",
"A4"
]
},
{
"cell_type": "code",
"execution_count": 89,
"id": "f737ddc7-6ada-4150-a67b-bc5042f29fa6",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\frac{\\left(1 - p\\right)^{n}}{p} + \\frac{1}{p}\\\\\\frac{\\left(1 - p\\right)^{n}}{p} - \\frac{4 \\left(1 - \\frac{p}{2}\\right)^{n}}{p} + \\frac{3}{p}\\\\- \\frac{\\left(1 - p\\right)^{n}}{2 p} + \\frac{8 \\left(1 - \\frac{p}{2}\\right)^{n}}{p} - \\frac{27 \\left(1 - \\frac{p}{3}\\right)^{n}}{2 p} + \\frac{6}{p}\\\\\\frac{\\left(1 - p\\right)^{n}}{6 p} - \\frac{8 \\left(1 - \\frac{p}{2}\\right)^{n}}{p} + \\frac{81 \\left(1 - \\frac{p}{3}\\right)^{n}}{2 p} - \\frac{128 \\left(1 - \\frac{p}{4}\\right)^{n}}{3 p} + \\frac{10}{p}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ -(1 - p)**n/p + 1/p],\n",
"[ (1 - p)**n/p - 4*(1 - p/2)**n/p + 3/p],\n",
"[ -(1 - p)**n/(2*p) + 8*(1 - p/2)**n/p - 27*(1 - p/3)**n/(2*p) + 6/p],\n",
"[(1 - p)**n/(6*p) - 8*(1 - p/2)**n/p + 81*(1 - p/3)**n/(2*p) - 128*(1 - p/4)**n/(3*p) + 10/p]])"
]
},
"execution_count": 89,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x0 = sy.Matrix([0,0,0,0])\n",
"b = sy.Matrix([1,1,1,1])\n",
"I4 = sy.eye(4)\n",
"S4 = (A4**n-I4)*(A4-I4).inv()*b\n",
"S4.applyfunc(lambda z: z.expand().powsimp())"
]
},
{
"cell_type": "code",
"execution_count": 90,
"id": "d69aab9b-f59c-438f-906a-9baf21d053da",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\frac{\\left(1 - p\\right)^{n} - 1}{p}\\\\\\frac{\\left(1 - p\\right)^{n} - 4 \\left(1 - \\frac{p}{2}\\right)^{n} + 3}{p}\\\\- \\frac{\\left(1 - p\\right)^{n} - 16 \\left(1 - \\frac{p}{2}\\right)^{n} + 27 \\left(1 - \\frac{p}{3}\\right)^{n} - 12}{2 p}\\\\\\frac{\\left(1 - p\\right)^{n} - 48 \\left(1 - \\frac{p}{2}\\right)^{n} + 243 \\left(1 - \\frac{p}{3}\\right)^{n} - 256 \\left(1 - \\frac{p}{4}\\right)^{n} + 60}{6 p}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ -((1 - p)**n - 1)/p],\n",
"[ ((1 - p)**n - 4*(1 - p/2)**n + 3)/p],\n",
"[ -((1 - p)**n - 16*(1 - p/2)**n + 27*(1 - p/3)**n - 12)/(2*p)],\n",
"[((1 - p)**n - 48*(1 - p/2)**n + 243*(1 - p/3)**n - 256*(1 - p/4)**n + 60)/(6*p)]])"
]
},
"execution_count": 90,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"S4.applyfunc(lambda z: z.expand().powsimp().factor())"
]
},
{
"cell_type": "code",
"execution_count": 91,
"id": "f7f7b765-23b8-408e-8a73-487e48b22edb",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\frac{\\left(1 - p\\right)^{n}}{p} + \\frac{1}{p}\\\\\\frac{\\left(1 - p\\right)^{n}}{p} - \\frac{4 \\left(1 - \\frac{p}{2}\\right)^{n}}{p} + \\frac{3}{p}\\\\- \\frac{\\left(1 - p\\right)^{n}}{2 p} + \\frac{8 \\left(1 - \\frac{p}{2}\\right)^{n}}{p} - \\frac{27 \\left(1 - \\frac{p}{3}\\right)^{n}}{2 p} + \\frac{6}{p}\\\\\\frac{\\left(1 - p\\right)^{n}}{6 p} - \\frac{8 \\left(1 - \\frac{p}{2}\\right)^{n}}{p} + \\frac{81 \\left(1 - \\frac{p}{3}\\right)^{n}}{2 p} - \\frac{128 \\left(1 - \\frac{p}{4}\\right)^{n}}{3 p} + \\frac{10}{p}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ -(1 - p)**n/p + 1/p],\n",
"[ (1 - p)**n/p - 4*(1 - p/2)**n/p + 3/p],\n",
"[ -(1 - p)**n/(2*p) + 8*(1 - p/2)**n/p - 27*(1 - p/3)**n/(2*p) + 6/p],\n",
"[(1 - p)**n/(6*p) - 8*(1 - p/2)**n/p + 81*(1 - p/3)**n/(2*p) - 128*(1 - p/4)**n/(3*p) + 10/p]])"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"S4.applyfunc(lambda z:z.factor().expand())"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "70805236-d935-452c-bedc-1f51574a0aa0",
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 92,
"id": "39a07040-85a8-46c0-8df3-744486a62420",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}\\left(1 - p\\right)^{n} & 0 & 0 & 0 & 0 & 0\\\\\\left(- \\frac{p - 2}{2}\\right)^{n} - \\left(1 - p\\right)^{n} & \\left(- \\frac{p - 2}{2}\\right)^{n} & 0 & 0 & 0 & 0\\\\\\frac{3 \\left(- \\frac{p - 3}{3}\\right)^{n}}{2} - 2 \\left(- \\frac{p - 2}{2}\\right)^{n} + \\frac{\\left(1 - p\\right)^{n}}{2} & 2 \\left(- \\frac{p - 3}{3}\\right)^{n} - 2 \\left(- \\frac{p - 2}{2}\\right)^{n} & \\left(- \\frac{p - 3}{3}\\right)^{n} & 0 & 0 & 0\\\\\\frac{8 \\left(- \\frac{p - 4}{4}\\right)^{n}}{3} - \\frac{9 \\left(- \\frac{p - 3}{3}\\right)^{n}}{2} + 2 \\left(- \\frac{p - 2}{2}\\right)^{n} - \\frac{\\left(1 - p\\right)^{n}}{6} & 4 \\left(- \\frac{p - 4}{4}\\right)^{n} - 6 \\left(- \\frac{p - 3}{3}\\right)^{n} + 2 \\left(- \\frac{p - 2}{2}\\right)^{n} & 3 \\left(- \\frac{p - 4}{4}\\right)^{n} - 3 \\left(- \\frac{p - 3}{3}\\right)^{n} & \\left(- \\frac{p - 4}{4}\\right)^{n} & 0 & 0\\\\\\frac{125 \\left(- \\frac{p - 5}{5}\\right)^{n}}{24} - \\frac{32 \\left(- \\frac{p - 4}{4}\\right)^{n}}{3} + \\frac{27 \\left(- \\frac{p - 3}{3}\\right)^{n}}{4} - \\frac{4 \\left(- \\frac{p - 2}{2}\\right)^{n}}{3} + \\frac{\\left(1 - p\\right)^{n}}{24} & \\frac{25 \\left(- \\frac{p - 5}{5}\\right)^{n}}{3} - 16 \\left(- \\frac{p - 4}{4}\\right)^{n} + 9 \\left(- \\frac{p - 3}{3}\\right)^{n} - \\frac{4 \\left(- \\frac{p - 2}{2}\\right)^{n}}{3} & \\frac{15 \\left(- \\frac{p - 5}{5}\\right)^{n}}{2} - 12 \\left(- \\frac{p - 4}{4}\\right)^{n} + \\frac{9 \\left(- \\frac{p - 3}{3}\\right)^{n}}{2} & 4 \\left(- \\frac{p - 5}{5}\\right)^{n} - 4 \\left(- \\frac{p - 4}{4}\\right)^{n} & \\left(- \\frac{p - 5}{5}\\right)^{n} & 0\\\\\\frac{54 \\left(- \\frac{p - 6}{6}\\right)^{n}}{5} - \\frac{625 \\left(- \\frac{p - 5}{5}\\right)^{n}}{24} + \\frac{64 \\left(- \\frac{p - 4}{4}\\right)^{n}}{3} - \\frac{27 \\left(- \\frac{p - 3}{3}\\right)^{n}}{4} + \\frac{2 \\left(- \\frac{p - 2}{2}\\right)^{n}}{3} - \\frac{\\left(1 - p\\right)^{n}}{120} & 18 \\left(- \\frac{p - 6}{6}\\right)^{n} - \\frac{125 \\left(- \\frac{p - 5}{5}\\right)^{n}}{3} + 32 \\left(- \\frac{p - 4}{4}\\right)^{n} - 9 \\left(- \\frac{p - 3}{3}\\right)^{n} + \\frac{2 \\left(- \\frac{p - 2}{2}\\right)^{n}}{3} & 18 \\left(- \\frac{p - 6}{6}\\right)^{n} - \\frac{75 \\left(- \\frac{p - 5}{5}\\right)^{n}}{2} + 24 \\left(- \\frac{p - 4}{4}\\right)^{n} - \\frac{9 \\left(- \\frac{p - 3}{3}\\right)^{n}}{2} & 12 \\left(- \\frac{p - 6}{6}\\right)^{n} - 20 \\left(- \\frac{p - 5}{5}\\right)^{n} + 8 \\left(- \\frac{p - 4}{4}\\right)^{n} & 5 \\left(- \\frac{p - 6}{6}\\right)^{n} - 5 \\left(- \\frac{p - 5}{5}\\right)^{n} & \\left(- \\frac{p - 6}{6}\\right)^{n}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ (1 - p)**n, 0, 0, 0, 0, 0],\n",
"[ (-(p - 2)/2)**n - (1 - p)**n, (-(p - 2)/2)**n, 0, 0, 0, 0],\n",
"[ 3*(-(p - 3)/3)**n/2 - 2*(-(p - 2)/2)**n + (1 - p)**n/2, 2*(-(p - 3)/3)**n - 2*(-(p - 2)/2)**n, (-(p - 3)/3)**n, 0, 0, 0],\n",
"[ 8*(-(p - 4)/4)**n/3 - 9*(-(p - 3)/3)**n/2 + 2*(-(p - 2)/2)**n - (1 - p)**n/6, 4*(-(p - 4)/4)**n - 6*(-(p - 3)/3)**n + 2*(-(p - 2)/2)**n, 3*(-(p - 4)/4)**n - 3*(-(p - 3)/3)**n, (-(p - 4)/4)**n, 0, 0],\n",
"[ 125*(-(p - 5)/5)**n/24 - 32*(-(p - 4)/4)**n/3 + 27*(-(p - 3)/3)**n/4 - 4*(-(p - 2)/2)**n/3 + (1 - p)**n/24, 25*(-(p - 5)/5)**n/3 - 16*(-(p - 4)/4)**n + 9*(-(p - 3)/3)**n - 4*(-(p - 2)/2)**n/3, 15*(-(p - 5)/5)**n/2 - 12*(-(p - 4)/4)**n + 9*(-(p - 3)/3)**n/2, 4*(-(p - 5)/5)**n - 4*(-(p - 4)/4)**n, (-(p - 5)/5)**n, 0],\n",
"[54*(-(p - 6)/6)**n/5 - 625*(-(p - 5)/5)**n/24 + 64*(-(p - 4)/4)**n/3 - 27*(-(p - 3)/3)**n/4 + 2*(-(p - 2)/2)**n/3 - (1 - p)**n/120, 18*(-(p - 6)/6)**n - 125*(-(p - 5)/5)**n/3 + 32*(-(p - 4)/4)**n - 9*(-(p - 3)/3)**n + 2*(-(p - 2)/2)**n/3, 18*(-(p - 6)/6)**n - 75*(-(p - 5)/5)**n/2 + 24*(-(p - 4)/4)**n - 9*(-(p - 3)/3)**n/2, 12*(-(p - 6)/6)**n - 20*(-(p - 5)/5)**n + 8*(-(p - 4)/4)**n, 5*(-(p - 6)/6)**n - 5*(-(p - 5)/5)**n, (-(p - 6)/6)**n]])"
]
},
"execution_count": 92,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Apn=Ap**n\n",
"Apn"
]
},
{
"cell_type": "code",
"execution_count": 93,
"id": "de9bc96f-8161-435c-89aa-b1c1d52cb71b",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\frac{\\left(1 - p\\right)^{n} - 1}{p}\\\\- \\frac{3 \\left(\\left(- \\frac{p - 2}{2}\\right)^{n} - 1\\right)}{p} - \\frac{\\left(- \\frac{p - 2}{2}\\right)^{n} - \\left(1 - p\\right)^{n}}{p}\\\\- \\frac{6 \\left(\\left(- \\frac{p - 3}{3}\\right)^{n} - 1\\right)}{p} - \\frac{3 \\cdot \\left(2 \\left(- \\frac{p - 3}{3}\\right)^{n} - 2 \\left(- \\frac{p - 2}{2}\\right)^{n}\\right)}{p} - \\frac{\\frac{3 \\left(- \\frac{p - 3}{3}\\right)^{n}}{2} - 2 \\left(- \\frac{p - 2}{2}\\right)^{n} + \\frac{\\left(1 - p\\right)^{n}}{2}}{p}\\\\- \\frac{10 \\left(\\left(- \\frac{p - 4}{4}\\right)^{n} - 1\\right)}{p} - \\frac{6 \\cdot \\left(3 \\left(- \\frac{p - 4}{4}\\right)^{n} - 3 \\left(- \\frac{p - 3}{3}\\right)^{n}\\right)}{p} - \\frac{3 \\cdot \\left(4 \\left(- \\frac{p - 4}{4}\\right)^{n} - 6 \\left(- \\frac{p - 3}{3}\\right)^{n} + 2 \\left(- \\frac{p - 2}{2}\\right)^{n}\\right)}{p} - \\frac{\\frac{8 \\left(- \\frac{p - 4}{4}\\right)^{n}}{3} - \\frac{9 \\left(- \\frac{p - 3}{3}\\right)^{n}}{2} + 2 \\left(- \\frac{p - 2}{2}\\right)^{n} - \\frac{\\left(1 - p\\right)^{n}}{6}}{p}\\\\- \\frac{15 \\left(\\left(- \\frac{p - 5}{5}\\right)^{n} - 1\\right)}{p} - \\frac{10 \\cdot \\left(4 \\left(- \\frac{p - 5}{5}\\right)^{n} - 4 \\left(- \\frac{p - 4}{4}\\right)^{n}\\right)}{p} - \\frac{6 \\cdot \\left(\\frac{15 \\left(- \\frac{p - 5}{5}\\right)^{n}}{2} - 12 \\left(- \\frac{p - 4}{4}\\right)^{n} + \\frac{9 \\left(- \\frac{p - 3}{3}\\right)^{n}}{2}\\right)}{p} - \\frac{3 \\cdot \\left(\\frac{25 \\left(- \\frac{p - 5}{5}\\right)^{n}}{3} - 16 \\left(- \\frac{p - 4}{4}\\right)^{n} + 9 \\left(- \\frac{p - 3}{3}\\right)^{n} - \\frac{4 \\left(- \\frac{p - 2}{2}\\right)^{n}}{3}\\right)}{p} - \\frac{\\frac{125 \\left(- \\frac{p - 5}{5}\\right)^{n}}{24} - \\frac{32 \\left(- \\frac{p - 4}{4}\\right)^{n}}{3} + \\frac{27 \\left(- \\frac{p - 3}{3}\\right)^{n}}{4} - \\frac{4 \\left(- \\frac{p - 2}{2}\\right)^{n}}{3} + \\frac{\\left(1 - p\\right)^{n}}{24}}{p}\\\\- \\frac{21 \\left(\\left(- \\frac{p - 6}{6}\\right)^{n} - 1\\right)}{p} - \\frac{15 \\cdot \\left(5 \\left(- \\frac{p - 6}{6}\\right)^{n} - 5 \\left(- \\frac{p - 5}{5}\\right)^{n}\\right)}{p} - \\frac{10 \\cdot \\left(12 \\left(- \\frac{p - 6}{6}\\right)^{n} - 20 \\left(- \\frac{p - 5}{5}\\right)^{n} + 8 \\left(- \\frac{p - 4}{4}\\right)^{n}\\right)}{p} - \\frac{6 \\cdot \\left(18 \\left(- \\frac{p - 6}{6}\\right)^{n} - \\frac{75 \\left(- \\frac{p - 5}{5}\\right)^{n}}{2} + 24 \\left(- \\frac{p - 4}{4}\\right)^{n} - \\frac{9 \\left(- \\frac{p - 3}{3}\\right)^{n}}{2}\\right)}{p} - \\frac{3 \\cdot \\left(18 \\left(- \\frac{p - 6}{6}\\right)^{n} - \\frac{125 \\left(- \\frac{p - 5}{5}\\right)^{n}}{3} + 32 \\left(- \\frac{p - 4}{4}\\right)^{n} - 9 \\left(- \\frac{p - 3}{3}\\right)^{n} + \\frac{2 \\left(- \\frac{p - 2}{2}\\right)^{n}}{3}\\right)}{p} - \\frac{\\frac{54 \\left(- \\frac{p - 6}{6}\\right)^{n}}{5} - \\frac{625 \\left(- \\frac{p - 5}{5}\\right)^{n}}{24} + \\frac{64 \\left(- \\frac{p - 4}{4}\\right)^{n}}{3} - \\frac{27 \\left(- \\frac{p - 3}{3}\\right)^{n}}{4} + \\frac{2 \\left(- \\frac{p - 2}{2}\\right)^{n}}{3} - \\frac{\\left(1 - p\\right)^{n}}{120}}{p}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ -((1 - p)**n - 1)/p],\n",
"[ -3*((-(p - 2)/2)**n - 1)/p - ((-(p - 2)/2)**n - (1 - p)**n)/p],\n",
"[ -6*((-(p - 3)/3)**n - 1)/p - 3*(2*(-(p - 3)/3)**n - 2*(-(p - 2)/2)**n)/p - (3*(-(p - 3)/3)**n/2 - 2*(-(p - 2)/2)**n + (1 - p)**n/2)/p],\n",
"[ -10*((-(p - 4)/4)**n - 1)/p - 6*(3*(-(p - 4)/4)**n - 3*(-(p - 3)/3)**n)/p - 3*(4*(-(p - 4)/4)**n - 6*(-(p - 3)/3)**n + 2*(-(p - 2)/2)**n)/p - (8*(-(p - 4)/4)**n/3 - 9*(-(p - 3)/3)**n/2 + 2*(-(p - 2)/2)**n - (1 - p)**n/6)/p],\n",
"[ -15*((-(p - 5)/5)**n - 1)/p - 10*(4*(-(p - 5)/5)**n - 4*(-(p - 4)/4)**n)/p - 6*(15*(-(p - 5)/5)**n/2 - 12*(-(p - 4)/4)**n + 9*(-(p - 3)/3)**n/2)/p - 3*(25*(-(p - 5)/5)**n/3 - 16*(-(p - 4)/4)**n + 9*(-(p - 3)/3)**n - 4*(-(p - 2)/2)**n/3)/p - (125*(-(p - 5)/5)**n/24 - 32*(-(p - 4)/4)**n/3 + 27*(-(p - 3)/3)**n/4 - 4*(-(p - 2)/2)**n/3 + (1 - p)**n/24)/p],\n",
"[-21*((-(p - 6)/6)**n - 1)/p - 15*(5*(-(p - 6)/6)**n - 5*(-(p - 5)/5)**n)/p - 10*(12*(-(p - 6)/6)**n - 20*(-(p - 5)/5)**n + 8*(-(p - 4)/4)**n)/p - 6*(18*(-(p - 6)/6)**n - 75*(-(p - 5)/5)**n/2 + 24*(-(p - 4)/4)**n - 9*(-(p - 3)/3)**n/2)/p - 3*(18*(-(p - 6)/6)**n - 125*(-(p - 5)/5)**n/3 + 32*(-(p - 4)/4)**n - 9*(-(p - 3)/3)**n + 2*(-(p - 2)/2)**n/3)/p - (54*(-(p - 6)/6)**n/5 - 625*(-(p - 5)/5)**n/24 + 64*(-(p - 4)/4)**n/3 - 27*(-(p - 3)/3)**n/4 + 2*(-(p - 2)/2)**n/3 - (1 - p)**n/120)/p]])"
]
},
"execution_count": 93,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Mpn=(Apn - sy.eye(6))*(Ap - sy.eye(6)).inv()*sy.Matrix([1,1,1,1,1,1])\n",
"Mpn"
]
},
{
"cell_type": "code",
"execution_count": 94,
"id": "c1df8e4b-ab55-4b8e-94d7-3e3111cb8267",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\frac{\\left(1 - p\\right)^{n}}{p} + \\frac{1}{p}\\\\\\frac{\\left(1 - p\\right)^{n}}{p} - \\frac{4 \\left(1 - \\frac{p}{2}\\right)^{n}}{p} + \\frac{3}{p}\\\\- \\frac{\\left(1 - p\\right)^{n}}{2 p} + \\frac{8 \\left(1 - \\frac{p}{2}\\right)^{n}}{p} - \\frac{27 \\left(1 - \\frac{p}{3}\\right)^{n}}{2 p} + \\frac{6}{p}\\\\\\frac{\\left(1 - p\\right)^{n}}{6 p} - \\frac{8 \\left(1 - \\frac{p}{2}\\right)^{n}}{p} + \\frac{81 \\left(1 - \\frac{p}{3}\\right)^{n}}{2 p} - \\frac{128 \\left(1 - \\frac{p}{4}\\right)^{n}}{3 p} + \\frac{10}{p}\\\\- \\frac{\\left(1 - p\\right)^{n}}{24 p} + \\frac{16 \\left(1 - \\frac{p}{2}\\right)^{n}}{3 p} - \\frac{243 \\left(1 - \\frac{p}{3}\\right)^{n}}{4 p} + \\frac{512 \\left(1 - \\frac{p}{4}\\right)^{n}}{3 p} - \\frac{3125 \\left(1 - \\frac{p}{5}\\right)^{n}}{24 p} + \\frac{15}{p}\\\\\\frac{\\left(1 - p\\right)^{n}}{120 p} - \\frac{8 \\left(1 - \\frac{p}{2}\\right)^{n}}{3 p} + \\frac{243 \\left(1 - \\frac{p}{3}\\right)^{n}}{4 p} - \\frac{1024 \\left(1 - \\frac{p}{4}\\right)^{n}}{3 p} + \\frac{15625 \\left(1 - \\frac{p}{5}\\right)^{n}}{24 p} - \\frac{1944 \\left(1 - \\frac{p}{6}\\right)^{n}}{5 p} + \\frac{21}{p}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ -(1 - p)**n/p + 1/p],\n",
"[ (1 - p)**n/p - 4*(1 - p/2)**n/p + 3/p],\n",
"[ -(1 - p)**n/(2*p) + 8*(1 - p/2)**n/p - 27*(1 - p/3)**n/(2*p) + 6/p],\n",
"[ (1 - p)**n/(6*p) - 8*(1 - p/2)**n/p + 81*(1 - p/3)**n/(2*p) - 128*(1 - p/4)**n/(3*p) + 10/p],\n",
"[ -(1 - p)**n/(24*p) + 16*(1 - p/2)**n/(3*p) - 243*(1 - p/3)**n/(4*p) + 512*(1 - p/4)**n/(3*p) - 3125*(1 - p/5)**n/(24*p) + 15/p],\n",
"[(1 - p)**n/(120*p) - 8*(1 - p/2)**n/(3*p) + 243*(1 - p/3)**n/(4*p) - 1024*(1 - p/4)**n/(3*p) + 15625*(1 - p/5)**n/(24*p) - 1944*(1 - p/6)**n/(5*p) + 21/p]])"
]
},
"execution_count": 94,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"MPn.applyfunc(lambda z: z.expand().powsimp())"
]
},
{
"cell_type": "code",
"execution_count": 169,
"id": "b7c12e1c-2554-4995-ba71-506d70f31f16",
"metadata": {},
"outputs": [],
"source": [
"abp = {p*a0:2*b1-(2-p)*b0-2,p*a1:2*b2-(2-p)*b1-2}"
]
},
{
"cell_type": "code",
"execution_count": 173,
"id": "e857581e-7170-49d0-965f-bb7f8384e407",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle b_{0} \\left(p^{2} - 3 p + 2\\right) + b_{1} \\cdot \\left(3 p - 4\\right) + 2 b_{2} - 3 p$"
],
"text/plain": [
"b0*(p**2 - 3*p + 2) + b1*(3*p - 4) + 2*b2 - 3*p"
]
},
"execution_count": 173,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ebp = (p*a1-(1-p)*p*a0-p).subs(abp).expand().collect([b2,b1,b0])\n",
"ebp"
]
},
{
"cell_type": "code",
"execution_count": 176,
"id": "5bd7dbc3-75d4-4695-8c92-d11d817cb33f",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle c_{0} \\left(p^{3} - 6 p^{2} + 11 p - 6\\right) + c_{1} \\cdot \\left(6 p^{2} - 22 p + 18\\right) + c_{2} \\cdot \\left(11 p - 18\\right) + 6 c_{3} - 6 p^{2}$"
],
"text/plain": [
"c0*(p**3 - 6*p**2 + 11*p - 6) + c1*(6*p**2 - 22*p + 18) + c2*(11*p - 18) + 6*c3 - 6*p**2"
]
},
"execution_count": 176,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bcp = {p*b0:3*c1-(3-p)*c0-3,p*b1:3*c2-(3-p)*c1-3,p*b2:3*c3-(3-p)*c2-3}\n",
"ecp = (ebp*p).expand().subs(bcp).expand().collect([c3,c2,c1,c0])\n",
"ecp"
]
},
{
"cell_type": "code",
"execution_count": 177,
"id": "68965679-7451-45e8-b8eb-03c45dc870ad",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle d_{0} \\left(p^{4} - 10 p^{3} + 35 p^{2} - 50 p + 24\\right) + d_{1} \\cdot \\left(10 p^{3} - 70 p^{2} + 150 p - 96\\right) + d_{2} \\cdot \\left(35 p^{2} - 150 p + 144\\right) + d_{3} \\cdot \\left(50 p - 96\\right) + 24 d_{4} - 10 p^{3}$"
],
"text/plain": [
"d0*(p**4 - 10*p**3 + 35*p**2 - 50*p + 24) + d1*(10*p**3 - 70*p**2 + 150*p - 96) + d2*(35*p**2 - 150*p + 144) + d3*(50*p - 96) + 24*d4 - 10*p**3"
]
},
"execution_count": 177,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cdp = {p*c0:4*d1-(4-p)*d0-4, p*c1:4*d2-(4-p)*d1-4, p*c2:4*d3-(4-p)*d2-4, p*c3:4*d4-(4-p)*d3-4}\n",
"dcp = (ecp*p).expand().subs(cdp).expand().collect([d4,d3,d2,d1,d0])\n",
"dcp"
]
},
{
"cell_type": "code",
"execution_count": 208,
"id": "56a4b3f6-c4d9-4166-8817-82eb782933f7",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0\\\\1 & p - 1 & 0 & 0 & 0 & 0 & 0\\\\2 & 3 p - 4 & p^{2} - 3 p + 2 & 0 & 0 & 0 & 0\\\\6 & 11 p - 18 & 6 p^{2} - 22 p + 18 & p^{3} - 6 p^{2} + 11 p - 6 & 0 & 0 & 0\\\\24 & 50 p - 96 & 35 p^{2} - 150 p + 144 & 10 p^{3} - 70 p^{2} + 150 p - 96 & p^{4} - 10 p^{3} + 35 p^{2} - 50 p + 24 & 0 & 0\\\\120 & 274 p - 600 & 225 p^{2} - 1096 p + 1200 & 85 p^{3} - 675 p^{2} + 1644 p - 1200 & 15 p^{4} - 170 p^{3} + 675 p^{2} - 1096 p + 600 & p^{5} - 15 p^{4} + 85 p^{3} - 225 p^{2} + 274 p - 120 & 0\\\\720 & 1764 p - 4320 & 1624 p^{2} - 8820 p + 10800 & 735 p^{3} - 6496 p^{2} + 17640 p - 14400 & 175 p^{4} - 2205 p^{3} + 9744 p^{2} - 17640 p + 10800 & 21 p^{5} - 350 p^{4} + 2205 p^{3} - 6496 p^{2} + 8820 p - 4320 & p^{6} - 21 p^{5} + 175 p^{4} - 735 p^{3} + 1624 p^{2} - 1764 p + 720\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ 1, 0, 0, 0, 0, 0, 0],\n",
"[ 1, p - 1, 0, 0, 0, 0, 0],\n",
"[ 2, 3*p - 4, p**2 - 3*p + 2, 0, 0, 0, 0],\n",
"[ 6, 11*p - 18, 6*p**2 - 22*p + 18, p**3 - 6*p**2 + 11*p - 6, 0, 0, 0],\n",
"[ 24, 50*p - 96, 35*p**2 - 150*p + 144, 10*p**3 - 70*p**2 + 150*p - 96, p**4 - 10*p**3 + 35*p**2 - 50*p + 24, 0, 0],\n",
"[120, 274*p - 600, 225*p**2 - 1096*p + 1200, 85*p**3 - 675*p**2 + 1644*p - 1200, 15*p**4 - 170*p**3 + 675*p**2 - 1096*p + 600, p**5 - 15*p**4 + 85*p**3 - 225*p**2 + 274*p - 120, 0],\n",
"[720, 1764*p - 4320, 1624*p**2 - 8820*p + 10800, 735*p**3 - 6496*p**2 + 17640*p - 14400, 175*p**4 - 2205*p**3 + 9744*p**2 - 17640*p + 10800, 21*p**5 - 350*p**4 + 2205*p**3 - 6496*p**2 + 8820*p - 4320, p**6 - 21*p**5 + 175*p**4 - 735*p**3 + 1624*p**2 - 1764*p + 720]])"
]
},
"execution_count": 208,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"k = 7\n",
"Tp = sy.zeros(k,k)\n",
"Tp[0,0]=1\n",
"for i in range(1,k):\n",
" Tp[i,0] = i*Tp[i-1,0]\n",
" for j in range(1,k):\n",
" Tp[i,j] = (i*Tp[i-1,j]-(i-p)*Tp[i-1,j-1]).expand()\n",
"Tp"
]
},
{
"cell_type": "code",
"execution_count": 189,
"id": "b3f44c70-3989-4c0c-9fce-b1bbc2158531",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}1\\\\p\\\\p^{2}\\\\p^{3}\\\\p^{4}\\\\p^{5}\\\\p^{6}\\end{matrix}\\right]$"
],
"text/plain": [
"Matrix([\n",
"[ 1],\n",
"[ p],\n",
"[p**2],\n",
"[p**3],\n",
"[p**4],\n",
"[p**5],\n",
"[p**6]])"
]
},
"execution_count": 189,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Tp*sy.ones(k,1)"
]
},
{
"cell_type": "code",
"execution_count": 206,
"id": "de147be4-8bb0-4d30-9320-928b75025d56",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle p + x - 1$"
],
"text/plain": [
"p + x - 1"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle p + x - 1$"
],
"text/plain": [
"p + x - 1"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[1 - p]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle p^{2} - 3 p + 2 x^{2} + x \\left(3 p - 4\\right) + 2$"
],
"text/plain": [
"p**2 - 3*p + 2*x**2 + x*(3*p - 4) + 2"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\left(p + x - 1\\right) \\left(p + 2 x - 2\\right)$"
],
"text/plain": [
"(p + x - 1)*(p + 2*x - 2)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[1 - p, 1 - p/2]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle p^{3} - 6 p^{2} + 11 p + 6 x^{3} + x^{2} \\cdot \\left(11 p - 18\\right) + x \\left(6 p^{2} - 22 p + 18\\right) - 6$"
],
"text/plain": [
"p**3 - 6*p**2 + 11*p + 6*x**3 + x**2*(11*p - 18) + x*(6*p**2 - 22*p + 18) - 6"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\left(p + x - 1\\right) \\left(p + 2 x - 2\\right) \\left(p + 3 x - 3\\right)$"
],
"text/plain": [
"(p + x - 1)*(p + 2*x - 2)*(p + 3*x - 3)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[1 - p, 1 - p/2, 1 - p/3]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle p^{4} - 10 p^{3} + 35 p^{2} - 50 p + 24 x^{4} + x^{3} \\cdot \\left(50 p - 96\\right) + x^{2} \\cdot \\left(35 p^{2} - 150 p + 144\\right) + x \\left(10 p^{3} - 70 p^{2} + 150 p - 96\\right) + 24$"
],
"text/plain": [
"p**4 - 10*p**3 + 35*p**2 - 50*p + 24*x**4 + x**3*(50*p - 96) + x**2*(35*p**2 - 150*p + 144) + x*(10*p**3 - 70*p**2 + 150*p - 96) + 24"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\left(p + x - 1\\right) \\left(p + 2 x - 2\\right) \\left(p + 3 x - 3\\right) \\left(p + 4 x - 4\\right)$"
],
"text/plain": [
"(p + x - 1)*(p + 2*x - 2)*(p + 3*x - 3)*(p + 4*x - 4)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[1 - p, 1 - p/2, 1 - p/3, 1 - p/4]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle p^{5} - 15 p^{4} + 85 p^{3} - 225 p^{2} + 274 p + 120 x^{5} + x^{4} \\cdot \\left(274 p - 600\\right) + x^{3} \\cdot \\left(225 p^{2} - 1096 p + 1200\\right) + x^{2} \\cdot \\left(85 p^{3} - 675 p^{2} + 1644 p - 1200\\right) + x \\left(15 p^{4} - 170 p^{3} + 675 p^{2} - 1096 p + 600\\right) - 120$"
],
"text/plain": [
"p**5 - 15*p**4 + 85*p**3 - 225*p**2 + 274*p + 120*x**5 + x**4*(274*p - 600) + x**3*(225*p**2 - 1096*p + 1200) + x**2*(85*p**3 - 675*p**2 + 1644*p - 1200) + x*(15*p**4 - 170*p**3 + 675*p**2 - 1096*p + 600) - 120"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\left(p + x - 1\\right) \\left(p + 2 x - 2\\right) \\left(p + 3 x - 3\\right) \\left(p + 4 x - 4\\right) \\left(p + 5 x - 5\\right)$"
],
"text/plain": [
"(p + x - 1)*(p + 2*x - 2)*(p + 3*x - 3)*(p + 4*x - 4)*(p + 5*x - 5)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[1 - p, 1 - p/2, 1 - p/3, 1 - p/4, 1 - p/5]"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle p^{6} - 21 p^{5} + 175 p^{4} - 735 p^{3} + 1624 p^{2} - 1764 p + 720 x^{6} + x^{5} \\cdot \\left(1764 p - 4320\\right) + x^{4} \\cdot \\left(1624 p^{2} - 8820 p + 10800\\right) + x^{3} \\cdot \\left(735 p^{3} - 6496 p^{2} + 17640 p - 14400\\right) + x^{2} \\cdot \\left(175 p^{4} - 2205 p^{3} + 9744 p^{2} - 17640 p + 10800\\right) + x \\left(21 p^{5} - 350 p^{4} + 2205 p^{3} - 6496 p^{2} + 8820 p - 4320\\right) + 720$"
],
"text/plain": [
"p**6 - 21*p**5 + 175*p**4 - 735*p**3 + 1624*p**2 - 1764*p + 720*x**6 + x**5*(1764*p - 4320) + x**4*(1624*p**2 - 8820*p + 10800) + x**3*(735*p**3 - 6496*p**2 + 17640*p - 14400) + x**2*(175*p**4 - 2205*p**3 + 9744*p**2 - 17640*p + 10800) + x*(21*p**5 - 350*p**4 + 2205*p**3 - 6496*p**2 + 8820*p - 4320) + 720"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$\\displaystyle \\left(p + x - 1\\right) \\left(p + 2 x - 2\\right) \\left(p + 3 x - 3\\right) \\left(p + 4 x - 4\\right) \\left(p + 5 x - 5\\right) \\left(p + 6 x - 6\\right)$"
],
"text/plain": [
"(p + x - 1)*(p + 2*x - 2)*(p + 3*x - 3)*(p + 4*x - 4)*(p + 5*x - 5)*(p + 6*x - 6)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"[1 - p, 1 - p/2, 1 - p/3, 1 - p/4, 1 - p/5, 1 - p/6]"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"x = sy.symbols('x')\n",
"for i in range(1,k):\n",
" px = sum(Tp[i,j]*x**(i-j) for j in range(i+1))\n",
" display(px)\n",
" display(px.factor())\n",
" display(sy.solve(px,x))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bd28dae0-4b75-4127-aea7-0acac96fadc2",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment