API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Sep 1st, 2015
70
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 7.70 KB | None | 0 0
raw download clone embed print report
  1. {
  2. "cells": [
  3. {
  4. "cell_type": "code",
  5. "execution_count": null,
  6. "metadata": {
  7. "collapsed": true
  8. },
  9. "outputs": [],
  10. "source": [
  11. "from numpy import *\n",
  12. "from numpy.linalg import inv, norm, det\n",
  13. "from scipy.linalg import sqrtm"
  14. ]
  15. },
  16. {
  17. "cell_type": "markdown",
  18. "metadata": {},
  19. "source": [
  20. "# The Polar Decomposition\n",
  21. "\n",
  22. "The polar decomposition, derivable from the singular value decomposition, states that the deformation gradient $F_{iJ}$ can be factorized as\n",
  23. "\n",
  24. "$$ F_{iJ} = R_{iK}U_{KJ} $$\n",
  25. "\n",
  26. "where $R_{iK}$ is unitary and $U_{KJ}$ symmetric-positive-definite.\n",
  27. "\n",
  28. "Left multiplying both sides by $F_{iL}$ gives\n",
  29. "\n",
  30. "$$ \\begin{align}\n",
  31. " F_{iL} F_{iJ} &= F_{iL} R_{iK}U_{KJ} \\\\\n",
  32. " &= R_{iM}U_{ML} R_{iK}U_{KJ} \\\\\n",
  33. " &= U_{ML} R_{iM}R_{iK} U_{KJ} \\\\\n",
  34. " &= U_{ML} \\delta_{MK} U_{KJ} \\\\\n",
  35. " &= U_{ML} U_{MJ} \\\\\n",
  36. " &= U_{LM} U_{MJ}\n",
  37. "\\end{align} $$\n",
  38. "\n",
  39. "In direct the notation:\n",
  40. "\n",
  41. "$$ \\boldsymbol{F}^{\\mathrm{T}}{\\cdot}\\boldsymbol{F} =\n",
  42. " \\boldsymbol{U}{\\cdot}\\boldsymbol{U} $$\n",
  43. "\n",
  44. "Analytically, the polar decomposition can be found from\n",
  45. "\n",
  46. "- Step 1\n",
  47. "\n",
  48. " $$\\boldsymbol{U} = \\left[\n",
  49. " \\boldsymbol{F}^{\\mathrm{T}}{\\cdot}\\boldsymbol{F}\n",
  50. " \\right]^{1/2}$$\n",
  51. " \n",
  52. "- Step 2\n",
  53. "\n",
  54. " $$\\boldsymbol{R} = \\boldsymbol{F}{\\cdot}\\boldsymbol{U}^{-1}$$\n",
  55. " \n",
  56. "A few red flags should be raised when viewing the preceding algorithm:\n",
  57. "\n",
  58. "- In Step 1, the matrix square root is taken. This is an expensive calculation in its own right. Often requiring an eigen analysis.\n",
  59. "\n",
  60. "- In Step 2, the matrix inverse is taken. For symmetric-positive-definite matrices (like $\\boldsymbol{U}$), there exist efficient procedures to determine the inverse. But, it's best to use algorithms that do not require it."
  61. ]
  62. },
  63. {
  64. "cell_type": "markdown",
  65. "metadata": {},
  66. "source": [
  67. "## Iterative Polar Decomposition\n",
  68. "\n",
  69. "The Taylor series expansion of $\\left(\\boldsymbol{I}-\\boldsymbol{X}\\right)^n$, where $\\boldsymbol{X}$ is a symmetric second-order tensor, is\n",
  70. "\n",
  71. "$$\n",
  72. "\\left(\\boldsymbol{I}-\\boldsymbol{X}\\right)^n \\approx\n",
  73. " \\boldsymbol{I} - n\\boldsymbol{X} + \n",
  74. " \\mathcal{O}\\left(\\boldsymbol{X}^2\\right)\n",
  75. "$$\n",
  76. "\n",
  77. "Let $\\boldsymbol{X} = \\boldsymbol{I} - \\boldsymbol{A}$, then\n",
  78. "\n",
  79. "$$\n",
  80. "\\boldsymbol{A}^n \\approx\n",
  81. " \\boldsymbol{I} - n\\left(\\boldsymbol{I} - \\boldsymbol{X}\\right) + \n",
  82. " \\mathcal{O}\\left(\\left(\\boldsymbol{I} - \\boldsymbol{X}\\right)^2\\right)\n",
  83. "$$\n",
  84. "\n",
  85. "For $n=-1/2$,\n",
  86. "\n",
  87. "$$\n",
  88. "\\boldsymbol{A}^{-1/2} \\approx\n",
  89. " \\boldsymbol{I} + \\frac{1}{2}\\left(\n",
  90. " \\boldsymbol{I} - \\boldsymbol{X}\n",
  91. " \\right) + \n",
  92. " \\mathcal{O}\\left(\\left(\\boldsymbol{I} - \\boldsymbol{X}\\right)^2\\right)\n",
  93. "$$\n",
  94. "\n",
  95. "Recall that \n",
  96. "\n",
  97. "$$\\boldsymbol{F} = \\boldsymbol{R}{\\cdot}\\boldsymbol{U}$$\n",
  98. "\n",
  99. "Let \n",
  100. "\n",
  101. "$$\n",
  102. " \\boldsymbol{C} = \n",
  103. " \\boldsymbol{F}^T{\\cdot}\\boldsymbol{F} = \n",
  104. " \\boldsymbol{U}^2 \\Rightarrow\n",
  105. " \\boldsymbol{U}^{-1} =\n",
  106. " \\left(\\boldsymbol{F}^T{\\cdot}\\boldsymbol{F}\\right)^{-1/2}\n",
  107. "$$\n",
  108. "\n",
  109. "So,\n",
  110. "\n",
  111. "$$\n",
  112. " \\boldsymbol{R} = \n",
  113. " \\boldsymbol{F}{\\cdot}\\boldsymbol{U}^{-1} = \n",
  114. " \\boldsymbol{F}{\\cdot}\\boldsymbol{C}^{-1/2}\n",
  115. "$$\n",
  116. "\n",
  117. "If the eigenvalues of $\\boldsymbol{U}$ are close to $1$, the Taylor series expansion is \n",
  118. "\n",
  119. "$$\n",
  120. " \\boldsymbol{C}^{-1/2} = \n",
  121. " \\boldsymbol{I} + \\frac{1}{2}\\left(\\boldsymbol{I} - \\boldsymbol{C}\\right) =\n",
  122. " \\boldsymbol{I} + \\frac{1}{2}\\left(\\boldsymbol{I} - \\boldsymbol{F}^T{\\cdot}\\boldsymbol{F}\\right) =\n",
  123. " \\frac{3}{2}\\boldsymbol{I} - \\frac{1}{2}\\boldsymbol{F}^T{\\cdot}\\boldsymbol{F}\n",
  124. "$$\n",
  125. "\n",
  126. "So,\n",
  127. "\n",
  128. "$$\n",
  129. "\\boldsymbol{R} = \\frac{1}{2}\\boldsymbol{F}\\left(\n",
  130. " 3\\boldsymbol{I} - \\boldsymbol{F}^T{\\cdot}\\boldsymbol{F}\n",
  131. " \\right)\n",
  132. "$$\n",
  133. "\n",
  134. "Using fixed point iterations, $\\boldsymbol{R}$ can then be refined by\n",
  135. "\n",
  136. "$$\n",
  137. "\\boldsymbol{R}_{n+1} = \\frac{1}{2}\\boldsymbol{R}_n\\left(\n",
  138. " 3\\boldsymbol{I} - \\boldsymbol{R}^T_n{\\cdot}\\boldsymbol{R}_n\n",
  139. " \\right)\n",
  140. "$$\n",
  141. "\n",
  142. "$\\boldsymbol{R}_0$ can be chosen as $\\boldsymbol{F}$."
  143. ]
  144. },
  145. {
  146. "cell_type": "markdown",
  147. "metadata": {},
  148. "source": [
  149. "## Function Definitions"
  150. ]
  151. },
  152. {
  153. "cell_type": "markdown",
  154. "metadata": {},
  155. "source": [
  156. "### Analytic Polar Decomposition"
  157. ]
  158. },
  159. {
  160. "cell_type": "code",
  161. "execution_count": null,
  162. "metadata": {
  163. "collapsed": false
  164. },
  165. "outputs": [],
  166. "source": [
  167. "def AnalyticPolarDecomposition(F):\n",
  168. " U = sqrtm(dot(F.T, F))\n",
  169. " R = dot(F, inv(U))\n",
  170. " return R, U"
  171. ]
  172. },
  173. {
  174. "cell_type": "markdown",
  175. "metadata": {},
  176. "source": [
  177. "### Numeric Polar Decomposition"
  178. ]
  179. },
  180. {
  181. "cell_type": "code",
  182. "execution_count": null,
  183. "metadata": {
  184. "collapsed": false
  185. },
  186. "outputs": [],
  187. "source": [
  188. "def NumericPolarDecomposition(F, niter=20, tol=1e-8, disp=0):\n",
  189. " I = eye(3)\n",
  190. " R = F\n",
  191. " absmax = lambda x: amax(abs(x))\n",
  192. " for i in range(niter):\n",
  193. " R = .5 * dot(R, 3 * I - dot(R.T, R))\n",
  194. " if absmax(dot(R.T, R) - I) < tol:\n",
  195. " break\n",
  196. " else:\n",
  197. " raise RuntimeError('Fast polar decompositon failed to converge')\n",
  198. " U = dot(R.T, F)\n",
  199. " if disp:\n",
  200. " print 'Fast polar decomposition converged ',\n",
  201. " print 'in {0} iterations'.format(i)\n",
  202. " return R, U"
  203. ]
  204. },
  205. {
  206. "cell_type": "markdown",
  207. "metadata": {},
  208. "source": [
  209. "## Evaluation"
  210. ]
  211. },
  212. {
  213. "cell_type": "code",
  214. "execution_count": null,
  215. "metadata": {
  216. "collapsed": false
  217. },
  218. "outputs": [],
  219. "source": [
  220. "# Define the deformation gradient\n",
  221. "F = array([[1, .5, 0], [0, 1, 0], [0, 0, 1.5]])\n",
  222. "assert det(F) > 0\n",
  223. "F"
  224. ]
  225. },
  226. {
  227. "cell_type": "code",
  228. "execution_count": null,
  229. "metadata": {
  230. "collapsed": false
  231. },
  232. "outputs": [],
  233. "source": [
  234. "# Verify\n",
  235. "Ra, Ua = AnalyticPolarDecomposition(F)\n",
  236. "Ri, Ui = NumericPolarDecomposition(F, tol=1e-5, disp=1)\n",
  237. "assert allclose(Ra, Ri)\n",
  238. "assert allclose(Ua, Ui)"
  239. ]
  240. },
  241. {
  242. "cell_type": "code",
  243. "execution_count": null,
  244. "metadata": {
  245. "collapsed": false
  246. },
  247. "outputs": [],
  248. "source": [
  249. "%%timeit\n",
  250. "R, U = AnalyticPolarDecomposition(F)"
  251. ]
  252. },
  253. {
  254. "cell_type": "code",
  255. "execution_count": null,
  256. "metadata": {
  257. "collapsed": false
  258. },
  259. "outputs": [],
  260. "source": [
  261. "%%timeit\n",
  262. "R, U = NumericPolarDecomposition(F, tol=1e-5)"
  263. ]
  264. }
  265. ],
  266. "metadata": {
  267. "kernelspec": {
  268. "display_name": "Python 2",
  269. "language": "python",
  270. "name": "python2"
  271. },
  272. "language_info": {
  273. "codemirror_mode": {
  274. "name": "ipython",
  275. "version": 2
  276. },
  277. "file_extension": ".py",
  278. "mimetype": "text/x-python",
  279. "name": "python",
  280. "nbconvert_exporter": "python",
  281. "pygments_lexer": "ipython2",
  282. "version": "2.7.10"
  283. }
  284. },
  285. "nbformat": 4,
  286. "nbformat_minor": 0
  287. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!