{ "metadata": { "name": "", "signature": "sha256:27607ffc81ddf6d6f18f9408ec86e0a7444d4bda8ea6ca4daea0ce0b45a94573" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "import pylab as pl\n", "import numpy as np" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "\u95ee\u9898" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "http://www.v2ex.com/t/148964\n", "\n", "\u73b0\u6709\u4e00\u4e2a m*n \u7684\u77e9\u9635\uff0c\u77e9\u9635\u5916\u56f4\u4e00\u5468\u7684\u70b9\u7684\u503c\u662f\u5df2\u77e5\u7684\uff0c\u77e9\u9635\u5f53\u4e2d\u6bcf\u4e2a\u70b9\u7684\u503c\u90fd\u662f\u4e0a\u4e0b\u5de6\u53f3\u56db\u4e2a\u70b9\u7684\u5e73\u5747\u503c\uff0c\u8bf7\u6c42\u51fa\u6bcf\u4e2a\u70b9\u7684\u503c\u3002" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "\u91c7\u7528\u8fed\u4ee3\u6cd5" ] }, { "cell_type": "code", "collapsed": false, "input": [ "M, N = 20, 40\n", "np.random.seed(42)\n", "m = np.random.rand(M, N)\n", "m[1:-1, 1:-1] = 0\n", "\n", "count = 0\n", "while True:\n", " count += 1\n", " avg = 0.25 * (m[:-2, 1:-1] + m[2:, 1:-1] + m[1:-1, :-2] + m[1:-1, 2:])\n", " err = avg - m[1:-1, 1:-1]\n", " if np.max(np.abs(err)) < 1e-15:\n", " break\n", " m[1:-1, 1:-1] = avg\n", "print count" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "3490\n" ] } ], "prompt_number": 11 }, { "cell_type": "code", "collapsed": false, "input": [ "pl.imshow(m, interpolation=\"nearest\");" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAWkAAADDCAYAAABAi+ByAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAE0tJREFUeJzt3XmQnMV5x/Hfo9WBjpUEEhKSURCRMQSHQrgIAhmsheA4\nRRBHjCExGIQpiIEyrviIi5CClUNwqjCJ7TLEVbgKH4C5BSYOxSFYkDiNAROwOYOIbBGEBCuttNLq\n6vyxh1arffvp2XlnpoW+n6qt2tHzTvcz/c48mp3p7tdCCAIA5GlYoxMAABSjSANAxijSAJAxijQA\nZIwiDQAZo0gDQMaGV3KwmY2WdLOkvSX9KITw4wFx5vMBwBCEEGywf7dK5kmb2VmS9pT0A0mLJf1F\nCGFLv3jYa8vv+47vXHiNxlzxtb7bVw7f1+3jwh/E469cvJ/bxp/svcw9ZgfrW6Wxrdtvr7rTvcum\nCadF42+2+3ketOjtaNz++ooB/9ImqaXv1hVa6PbROjset5kJ53+Zf8gO3mmVprX2S8K/y3UnL4jG\n79apbhsP3HRyNL7szCk73P731vX6+9axfbdn2kq3D90YD4875T23iflj743Gbz73vB1utz4vtR62\n/bYt88/ZNY9cFI3/uR5y25ht/xiNvxku3eH291o79JXW5r7bs/Z/x+1Dy+6Lhq/TCW4T7zvxy/o1\n0fqa1PqxHeN2UcJr4MS3ouGw9I/dJuZ/8rbC2H/a6YVFutKPO2ZLWhq6K/tLkg6o8P4AgApUWqTH\nS1rX8/s6SRPLTQcA0F+lRXqtpN6/Z8ZJWhM7eMS8o4aSU/2NaGl0BolmNjqBNONaGp1BkiNbRjQ6\nhSQt+zQ6gzRzWkY2OgVXy6RGZ1C5ir44lPQrSUeb2fOSDpX02sADOhde0/f7LlOkR7Y0OoNEMxud\nQJrmlkZnkOSoXaCoSFLLtEZnkObIllGNTsGVS5Fe1fayVre9nHRspUX6Lkk/k3S6pJ/2/9KwV/8v\nCgEAO5vc8nFNbvl43+3XF95ReGxFRTqEsFnS38SOeX/4lMLY3OB/z/gHrY7Gl+gYt41xy+LfsK97\nae94A/uc5PaxZL/4Xwmb5L9Le/XUA6PxeeHIaPyrW/0/159r+mg0vp9ecdt4e/FB8QM2xsPDDl/v\n9nFr/GmlDRrttrHvma9H4xfo+ngDz7pdSPvEH+y0sSsSGom7/oazovFD9ZTbxrM6vKq4JOmlBdHw\n+Zocv/9S54khaZ+PxF8DFyZMuNGj8fD6k+Kf6p4+6iduF79onx+Nz52w2G3jSYu/nouwmAUAMkaR\nBoCMUaQBIGMUaQDIGEUaADJGkQaAjFGkASBjFGkAyFhFW5W6jZmFbXtF4v/lt/E/c+IbFbTpWLeN\nh3R8NL5C06PxZnW4fczVE9H4RH3gtrFJ8WW072pqNP5RveH2sUrxdbCLnbGSpN/q4Gi8o6s5Gm8e\n5Y/nFL3rHuPxxnOkuqLxlc54d7exKRqfrFVuG9MUX/Ay2VnQ5T0vJKlJOy0G3oE3VpL/Opmo9mh8\nq5rcPmY6++CepuKVeL2814E3Xk9rjtvHCsXX5k+Vv+om9li/YHeWtlUpAKCOKNIAkDGKNABkjCIN\nABmjSANAxijSAJAxijQAZKzSK7O4ovtaR+ZQ99qgMdF4yhzmA/VqND4pYS6rx9vUv117Vt3HKGde\n73LNcNtod64VnDKeM7Q8Gu8YFZ8n3aStbh/DE47x2+is6v6TnPnJkj//eLQ2uG1sdV52HYqPp/e8\nkKQtzhzllHMy3ZnP7bUxJuF8eOsJUp7jnc4FIbzxTlnTMNWZxz9bL7htfOJHvyuMfSFyP95JA0DG\nKNIAkDGKNABkjCINABmjSANAxijSAJAxijQAZIwiDQAZK30xy8S73imM3TjqTP/+zkbiKZPwvUUJ\nXhtdzkKVFN5ilxRdzsbsKZuqe214m9hL/jnx2kjJ05Ny3qvtJ2UsvDxSFpp4yhgvb3FQGYuHvPFK\nWcwyymnDW4wl+ePlnTNvgVJZbehZ/5DB8E4aADJGkQaAjFGkASBjFGkAyBhFGgAyRpEGgIxRpAEg\nY6XPk35v47TC2IbhI9z7r2iaHo2vczZEl/yN7MuYw1zGXFZvY3avjzIeR8r8Y2++q9dGyliVMZ71\n4M+X9cfTm1Ob0obfR+3bqMfjSFmz0ORcKCRpDrNjtSZH4/frM24bP3/0i5GoFUYqfidtZgeZ2dtm\n9oiZXVXp/QEA6YbyTnq8pB+GEL5ddjIAgB0N5TPpCZLWlJ0IAGBnQ/3i8Dwze9zMzi81GwDADir+\nuCOE8KCkB81shKTHzOyWEELfN3Xf6vchyLyjpXnHlJEmAHyIrG+TOtuSDq24SJvZsBDCthDCZjPr\nlLStf/zySyttEQB2M2Nbun96rV5YeOhQvjj8rJldrO45I3eGENYPoQ0AQIKhfNxxu6Tba5ALAGCA\n0hez3BTZo3vBDZvd+3ctiG/Yv0LxxS5S9QsOythg3luo0t1G6cNfsZSJ/imb4cekjIUnh7GS/PFK\n2Uy/jOdftX2ktRF/rN6G/WW8jjY5F61I4Z2TlDxHOhdzOEO3um3cfMh5hTF7pfh+LAsHgIxRpAEg\nYxRpAMgYRRoAMkaRBoCMUaQBIGMUaQDIWOmTT+fFgvv791+tSdF4R8Km/xs0OhrPZYN5L4965Jky\nr1fOHFFPUynzpKuf91uG+mymX/s5zmXM5642hxQprwFvLnUZzx2v7tysw902Tpl1fyRa4qb/AID6\noUgDQMYo0gCQMYo0AGSMIg0AGaNIA0DGKNIAkDGKNABkzEII5TVmFm4LJxbGUybHtyty1QBJnRrj\nttHpLGbxJ7/7E+i9jexTNqmvdjFLSp5dJTzWai9wUMZY7CpyWYji91GPRTn+YpZcHovHu/DFRH3g\ntjFHzxTGWuwZhRAGXdHCO2kAyBhFGgAyRpEGgIxRpAEgYxRpAMgYRRoAMkaRBoCMlb7p/zd0dWFs\nlt5w77+/lkXjE9XutuHNadydlLHxerVS5sJ6c60/TKqdG1yPOc5pbTT+uSVJXRpZ8z68efz76y23\njXlXFM+TjuGdNABkjCINABmjSANAxijSAJAxijQAZIwiDQAZo0gDQMYo0gCQMXfTfzM7VtIiSXtJ\nOlLSdyRtlHR2COH3A44N0vLixtr2dRM6Yt6j0fgsvem24S14GaPOaLwem4inKGPT/zIWiXib9tdj\nw/5cLgpQj+dGffqofiFKGRv2e1Kev9U+P1MuJLJSU6LxJ9bPddtYN+5/I9HDh7bpv5k1Szpb0ouS\nTNJlkuZL+rqkf3CzAgBUJVqkQwgdIYRzpb7/MseEEFaHEJ6TdHDNswOA3Vyln0n3P35EmYkAAHZW\n6QZL2/r9XvCB1L/1+/2onh8AwHbPSvp10pGVFum1ZjZZ0gxJrw5+yFcrbBIAdjeH9/z0ur7wyNQi\nHXp+Fqp7psdWSQuGlhwAIFVSkQ4hHNfz63OSjqldOgCA/krf9F96tzi0hz9P2psXmTJftoz5xZ4c\n5lKXkUPKPNQcNnfPYbzrZVeZw7yr8F7vm0q4aMDssS+4x5wRbi2MfXnQGdLdWHEIABmjSANAxijS\nAJAxijQAZIwiDQAZo0gDQMYo0gCQMXc/6YoaMwtrthTvu/StpsvdNt7V1Gg8ZQ7pGG2Ixkepy+nD\nn2PqHZPSRrXztcuYM17GftP1UMa8393pseaw73UZ6xHK2DN9g7NfdIea3T5WaHo0/vTqOW4bm2eO\nLw6us6HtJw0AaCyKNABkjCINABmjSANAxijSAJAxijQAZIwiDQAZo0gDQMZK3/T/raaZhbEZWu7e\n35u83qVRbhvegpcyFqKUsVigHosBPCkLJ7zFAh+WxRcp6nFOynh+VttHGcpY0JX2WqzunKQscpqk\nVdH4CZN+6bbxw44LC2PT2PQfAHZNFGkAyBhFGgAyRpEGgIxRpAEgYxRpAMgYRRoAMlb6POlDf/F6\nYWzxSce79/fmQafMU612nvRI56IAUh6b0Ndjzm4Z6nURBe+8b3We7vXapL7ax1rGWOSiHnmO1KZo\nPOW13K6J0fg9fzjVbeOefddFohcURngnDQAZo0gDQMYo0gCQMYo0AGSMIg0AGaNIA0DGKNIAkDGK\nNABkzF3MYmbHSlokaZKkf5HUImmDpCtCCI/tdPxXQmFbe8x7303oUxOWRONTtNJtw1uM4i0GGOVM\nfk/pw1s40Z1HfOFDZ/lrjXZS7YKalDZSFl9445lyTrx+NmlkNJ5yQYkcNv0flbDYqto+UuQwFinH\njFFn1Xl4dWf6R95x27gjnFYYey+y6X+0CphZs6SzJb3Y80/jJZ0aQvAzAgBULfpxRwihI4RwrtT3\nX9UESWtqnhUAQFLln0lvkXSPmd1vZgfWIiEAwHaVfOgZQgjnSJKZHSbpSkmf2+moD1q3/75HizS6\npYr0AODDZ1Pbk9rc9lTSsRV9M2VmTSGErZLaJW0b9KA9WytpEgB2OyNbjtLIlqP6bncu/G7hsalF\nOkgySdeb2cye25cMOUMAQJKkIh1COK7n1y/WMBcAwADlT8Rd1loY2vh/xbFeXRP8uaoeb46yN890\notrdPrw2Uubcdmq0Ex/jtuEp5yIK8XmoXhspm6p786BTzonHmyedMme8S+Oi8ZT58d5zx4t7m9in\nSNls3ztvZTw/PSkX4PCeO9O1IhqfpNVuH5O0Kho/dEXxxU56XXvI1wpjsRkcrDgEgIxRpAEgYxRp\nAMgYRRoAMkaRBoCMUaQBIGMUaQDIGEUaADJW+mKWcNzCwtjtB/7avf8jOjYaT5mE3649o/EZWl5V\nXPI3El+uGW4b3mNZqanReMqihhWa7rThLxZoVkc0vlqTovGUBTOz9GY07i1ISNGh5qrikv/c8hbM\nSP54emORwlsoNUYb3Da8PD/QxGi8jAs1pPAW1Xiv54/pVbeP5jUb4wf4JUN2VST4peIQ76QBIGMU\naQDIGEUaADJGkQaAjFGkASBjFGkAyBhFGgAyVvo86RsfLo5NduZdSv6c2zl62m3jap1WVRtz9YTb\nx94r10Xj/zHlHLeNg/W7aPwO53Gcqx+7fdyqM6LxQ/Si24Y3Xnfr1Gh8it51+zhDt0bjs/SG28YG\nZ76sN5/Wm1MuSf+tQ6LxlDnh3hz7r+s70fhvdbDbx/Oa7eTgz5Nu0SPR+CKdEo2nvI48KReMuFcn\nReOXriy+fqAkaXFCIkc78ZRp/EsTjhkE76QBIGMUaQDIGEUaADJGkQaAjFGkASBjFGkAyBhFGgAy\nRpEGgIxZCKG8xsyCxkXauyOhkYfi4bC3+Xl883vR+OJwSzR+3DefdPtYG+9CE17e5rYRzon/H2mP\nr4rf/4DJbh/2evz8HhB+47bx2unxhRF2n/Mciq/JkSSFec55TVh2tXl+PD7CuebEz4872e3j85fe\nHT/gS5vdNnTiiGj4qZfiYzEs/KnbxTn6STR+RMKisE/aRdH4BRPj5/2ZD+ILfyTpIR0fjX9eN7lt\nzHTyXD+2+EIkknTfercLjQ/x1SyffsBfqXLfZ4pjJ0gKIQx64nknDQAZo0gDQMYo0gCQsdoW6S1t\nNW2+LC+0rWl0ComGuENLve0i573thUZnkMa/fHMmNrc1OgPXY9Vf97bualukt7bVtPmy/KZtbaNT\nSPR4oxNIs4uc9zb/O9MsPNfoBFLtAv85L6FIAwDKRJEGgIyVP08aAFCxonnSpRZpAEC5+LgDADJG\nkQaAjNWsSJvZaDNbZGZLzWxBrfqplpkdZGZvm9kjZnZVo/MZyMyONbN2MxtmZnPN7Akze9jM9m10\nbv31y7PJzP7VzJ7qGdNPNTq3XmY2x8yW9OR2Wq7jOUie2Y2nmf1lz2v7KTP7u4zHcmCe2Y2lK4RQ\nkx9JZ0n6siST9LCk4bXqq8o8j5B0aaPzKMitWdINkh6T1CTpl5ImSfqEpO83Or9IntdJmtbovAbJ\nc5qksZJGqnvSea7jOTDPa3MbT0lTes71MEnPZjyWA/PMbiy9n1p+3DFb0tLQPVIvSTqghn1VY4Kk\nLJcchhA6QgjnSn3XtR8TQlgdQnhO0sENTG0Hg+SZ5ZiGEN4JIfTuebZB0uhMx3NgnhOV2XiGEFaG\nELZKmi5ppfIdy4F5ZjeWnloW6fGS1vX8vk7dg5Or88zscTM7v9GJOPqfr/h+l421RdI9Zna/mR3Y\n6GQGcbGkn6n7HVavHMezN88sx9PMrlb3XgXXKuOxHJBnlmMZk7BL75CtVfefwZI0Tpn+7xVCeFDS\ng2Y2QtJjZnZLCKGj0XkV6L9Jda4LXEMI4RxJMrPDJF0p6XONTWk7M/szSUeEEP52wHclWY1n/zyl\n7s2hcxvPEMI3zOyfJD2g7o9nemU1lv3yvF/SX4UQ1uc2ljG1LNK/knS0mT0v6VBJr9WwryEzs2Eh\nhG0hhM1m1qkdC2Fu1prZZEkzJL3a6GSKmFlTz5+Y7cpoPM1sqqR/1vbLEGQ5ngPzzHE8zWyPEMLG\nEEJXzyK2NZmOZf88pe3v+LMZS08ti/Rd6v5T7XRJPw0hbKlhX9X4rJldrO4vOO/s91lgTkLPz0JJ\ni9T9TmVBIxMqENQ9jteb2cye25c0MqEBLpf0R5Lu7XnBXqY8x3Ngnm9lOJ6XmNl8dX8Ed5ukJcpz\nLAfm+d0MxzKKFYcAkDEWswBAxijSAJAxijQAZIwiDQAZo0gDQMYo0gCQMYo0AGSMIg0AGft/OSaA\nam5gC4kAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 12 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "\u91c7\u7528\u7a00\u758f\u77e9\u9635\u89e3\u7ebf\u6027\u65b9\u7a0b" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy import sparse\n", "from scipy.sparse import linalg" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 13 }, { "cell_type": "code", "collapsed": false, "input": [ "idx = np.arange(M * N).reshape(M, N)\n", "\n", "A = sparse.dok_matrix((M*N, M*N))\n", "b = sparse.dok_matrix((M*N, 1))\n", "i = 0\n", "for row in range(M):\n", " for col in range(N):\n", " if row == 0 or col == 0 or row == M-1 or col == N-1:\n", " A[i, idx[row, col]] = 1\n", " b[i, 0] = m[row, col]\n", " else:\n", " A[i, idx[row, col]] = 4\n", " A[i, idx[row-1, col]] = -1\n", " A[i, idx[row+1, col]] = -1\n", " A[i, idx[row, col-1]] = -1\n", " A[i, idx[row, col+1]] = -1\n", " i += 1\n", " \n", "\n", "x = linalg.spsolve(A, b)\n", "x.shape = M, N\n", "\n", "np.allclose(x, m)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 14, "text": [ "True" ] } ], "prompt_number": 14 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }