{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "020f3c56",
   "metadata": {},
   "source": [
    "线性方程组与高斯消元"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "db5e8325",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "解 x = [2. 1.]\n",
      "验算 A @ x = [5. 5.]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['font.sans-serif'] = ['SimHei']\n",
    "plt.rcParams['axes.unicode_minus'] = False\n",
    "\n",
    "# 1. 定义系数矩阵 A\n",
    "A = np.array([[2.0, 1.0],\n",
    "              [1.0, 3.0]])\n",
    "\n",
    "# 2. 定义常数向量 b\n",
    "b = np.array([5.0, 5.0])\n",
    "\n",
    "# 3. 使用 np.linalg.solve 求解线性方程组 Ax = b\n",
    "x = np.linalg.solve(A, b)\n",
    "print('解 x =', x)\n",
    "\n",
    "# 4. 验算：将求得的解 x 代回原方程 A @ x，看是否等于 b\n",
    "print('验算 A @ x =', A @ x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "981939a5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAHSCAYAAAAnuM2rAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAANeVJREFUeJztnQeYFEX6xr9dFpYclyDKETyCIEg8EBBFUZEgcgqCp38BMXCCnoQTzsihBD1A8QAVSSqewomAnAgeQTGAoBySw7IgkkTCsssmlu3/8xb2ODM7uzt5urven0/LzkxPT01Nv/VVffXVV3GGYRhCCNGG+FgXgBASXSh6QjSDoidEMyh6QjSDoidEMyh6QjSDoidEMyh6QjSDoieOBXFneXl5Pp/XGYq+CM6fPy8tWrSQjz76SI4fPy5dunRR/0aKFStWyOOPPy5ZWVkSbSCQQD8XAsrMzCzw9Q8//FDmz58vubm5IZfv/fffl0WLFvl9/urVq6Vx48ayfft2j+dvuOEGufnmm9XfW7du9dkwmHzyySfSunVrCRfr16+Xq666Sv3OMQNhuMQwdu7caWzZssV1nDx5Uj2/cuVKmAVj27ZtRnZ2tnHllVcagwcPjkgZ8vLyjPbt2xu9e/f2eP6nn34ySpQoYVSsWNGoUqWKxzFy5Miwff6+ffvUdw3myMzMzHe9jIwMo27dusaNN95opKSk+DxOnDiR732pqanGuXPnjJycHOPChQuuo23btka3bt08nsM558+fV+/xBvWIz87NzfV4/o477lCvpaenG7Vq1TLGjBlTYJ0sWrTIKFOmjOvxjz/+WGg9zJgxo9A6xvctXry40bNnTyNWJMSuubEWvXv3ll9++UXKli0rx44dkylTpsijjz4qS5culTZt2sgVV1whGRkZ8re//U1KlSolp06dkpycHHV+uXLlwlKG6dOny9dff63+jouLU/9OmzZN/vznP8uBAwfk008/leeee85lufr27SslSpSQcFGnTh05fPiwlClTxvX5r732mixbtkw+++wz9fi+++6TSpUqqXKBixcvqt5QYmJivus988wzkpKSoo66devmez0+Pl5ZvFtuucXj+f/7v/9T9V4QxYsXz/fc9ddfL+vWrXM9/uGHH1R9oa7S0tLUb1SsWDH1WoUKFeTChQvqe7799tvywQcfKGuP8niD9yQkJOT77JUrV0qDBg08zu3YsaNH2VJTU6V8+fKuugTVqlWThx56yOcQA8+hrLinfJUlXFD0v4If6x//+IcMGDBAdefwGDc0uqf48SAIcO7cOSldurT6G6KfMGGCjBgxIuTPh7CefPJJ1X296667ZNiwYeqz8C+4/PLL1Y122223ScWKFdVzZ86ckerVq0u4wM2Nxs0diAHiMD8TNzBuZPMxqFKlSr5rQUivvPKKuskfeeQRGTt2rOs1NALt27eX4cOH5xM8mDt3rhoOoCFxv/khsjFjxsjAgQM9hIJz3UWEvx977DEZOXKk1KtXT52PRuCbb75xNZJ4DxrS06dPq0b8mmuukXfffVf9C15//XVVH+j+43dGmdBQNG3aVL1es2ZN1z3hXn/uokcdeTca7syePdvjMT4HZYfRqVGjhkQKiv5XvC0mHmMcf+LECRk3bpw8/fTT6vmkpCT597//rcaF4QKtO6x5o0aNpE+fPh6vwRI9+OCD0rlzZ1m8eLF8++23smHDBiWo//3vf6rRCSf4vu5j9KNHj6qb9uDBg+oxGkDcmOZjU/TuvR2UE43njBkz1NgZfhDcyFOnTlXCQ28BDcGoUaN8lgE9CYAeBIRmNj4///yzqn93/wAaIfN8k9GjR6tyQ/R79+6VBQsWqF7UO++8I7t27ZK1a9fKjz/+qMb8LVu2lGbNmsnQoUM9GlA0vvhMvB9lwHvxuWYjXBBmbwKg14QGBc+99dZbqpx33nmnx/kvvfSS1K9fX/U04U9B/eI7RpSYDSwsRqtWrYy5c+d6/N2uXTsjLi7OGDdunOs8jKPXrl2rxokYd4aLs2fPGhMmTDB69eplHD58WI0PwQMPPGA88cQTxp///Gfj2WefVc/h9T59+hivv/668gOEk5tvvjng8fzMmTPzjYNnz57tegz/SJMmTYykpCQjPj7emD9/vl9lue2224r8bIyP3UF91KtXT72Gz8KB8TPG//hdhw4datx+++1Gx44d832e6cdxBz6WChUquB6fOXOm0PK88847+a6xePFiIyEhQfkk3H0f8Cnccsst6h7761//qsoYDSj6AkQPB1DVqlWNm266qdAfOS0tTTn44EzCv/6K8OLFi+oGgCPKBKLv0KGDsXDhQnVtNC5wXD399NPqul9++aU6r1OnTsaqVasiUg8QyIgRI1yPr7nmGuOVV15Rf6MM5cqVMz777DPX63Bszpkzx+e1UCfLli1TIvvd736nvh9ufIh/4MCBSvybN282Tp065VEP7g63xx9/vMCyfvzxx6o83mzYsMHYvXu3ceDAAdVIw2HozrRp01Qj9MMPPyiR4jPwuGTJkup3cS8/GhWIftOmTR6i//rrr9Xf7gecgt6if+ONN5Tge/Tooa7nzc8//2y8+OKLqnG69tprjYMHDxqRhlN2BXDttdfKG2+8oZwqTz31lBo/46hcubJ8/PHHyun3008/qfEuurH41xyDostZ1IEuH7p+6Hq6s2fPHuVEBM8//7x8//336twjR44oR9EXX3yhuo3mWNfXVNju3bvVYXaNA8F9DI2pSYyFGzZsqB7v27dPDUWaN2/u8R7TUYXuMLrwGI5gzI7u7MSJE+X2229XdYbhCIZG6DrD4KBe4T/B8ADj9fT0dI/r4nu/+uqrBdZhz549fTq82rZtq8r8xBNPqC6+OfZGfd5zzz0yefJk2bFjh9x4440yb948NXzBkGDTpk0eTjf4WVCHGO7gs9AVN8FwBmN298O7LPjshx9+WAYPHqyu369fP+Uncgd1g6EGHJf4N5LTwSYc0xfA73//ezXOwhxzyZIlPRxZaAjcnVdoCK6++mrlxMHhfuMUBH58OG7cHWLgpptuUg7FWrVqKW80bhjw3nvvqX9xo+K9pgOsW7duSlDuYB7YFCm+R7BgHAqnHcbg//nPf9R4HDz77LNK3N4ee9z0qCs0RBizYzx/2WWXqde+++475dybNGmS8ofgwPeAMw3lRAOBevUGMQvwX/hi+fLlcu+99+Z7Hg3KCy+8IPv371cNDGZdUNdoLFu0aKHe06NHD/nqq69UYwMxw1/QpEkTj9/uzTfflNq1a6vGfc6cOUqgGP8DNATejZS3Rx6NDO4L+DfQyKC8//rXv1xlhmMW/hn4i7p37y7Jyck+Z0HCTsT7Ejbt3ptdVoyx0cX905/+pI7ExESjS5cuRr9+/dRrvsaBwYL5Ynw2uvP4aV5++WXjD3/4g3qMcSrGhmaXGl1/dBd9fb459MC8e6DgO6F7f+zYMdWtff/9943hw4cbd955p6qHvn37Go0bNzZat26t5qxRFrPeTLp37x6QT+DRRx/1WRZ8JrrvtWvX9nlUq1bNY7xtMmvWLHVdvBdz7Ph74sSJHuc0atRIlXvJkiVGsWLFjP79+6suusk333yjxtp4H14H8KugPgId05vcfffdqtwYy2NIgxiGAQMGGNGGlt4NTO2Y00HuUVqwAJhGM7t8nTp1UtYUrb05fRcOMH0EjzW61GDz5s1qGAHrA2sOL6876JZilgGW1X1aKNQwU7wfswiwiogFgMcZFmrQoEGqGwrrD488PNq+gIVEtxhW3/RmY74cMxD4bqb1xznZ2dlqaFQQ+MxALT16GOgtwcpiyg7Ta96f0b17dzWNCCuOKVf0QExQLkz5oUtvDm0AZgDOnj3riqzDbIs7mAkoDHwP9BQwU4NuPMqEOIioE/VmxqI0b97cmD59umrt8bcZWQXL58t7HwngPJw6dWo+7z0O0xsNq2P+bVqX9evXe1xn165d6vDlHPPX0n/66afG0aNH1XPHjx9XVvUvf/mL6zxYK+DL0pts3LhROePgsESUI8qK77Zjxw71GVlZWYWWBZa+qF6CL0tvAm/49u3b1UwCoihNS56Xl2cMGzZMvR89NhPTCYtywYEKB+BHH33ksvTujjxEaHoDK16YpQdvvvmmej+ce/iNYgEdeb+C8SWsNsbYGC97B174soYYJ3qP64IF1nrLli35rLkZmIM5XJQR51155ZVqjhmP8flwOrqDXgiOQ4cOBVwOs5dw6623uiwy5q8RSIJ4BRPTchbkLIRTDIFEmJ/HmNqdnTt3quCXVq1aybZt2wotD3wDphPV+0Asvi8Q9ANfBhyl7dq1U70UBLtgvHz48GFVLgTioPeGqD0EC6HHgZ4B/BY4Dz0rX/dAYXH6oKA1BqgnOIbh3IPzEo7A/v37q8+P9gIgdu9/BT+6iXtXDz8ywklxmKCb6o4ZOhkKcOogMAOeZ9wEHTp0UF3pzz//XHVh3R085k0Cx1lhXeNgKOimhePLBF54BLggAAYiqlq1qse58IijC4v3wBEK8cFhZwKxQfDwZkOUmP24//771WsnT55U18WwBfWK71eQRxvCx++DmQqUG0KC4w0NFp7v1auXXHfddWrog8+H4ObMmaO62JgVgbMUZURdozyYOUEEIXCfHXH3uJt/m5F5hdUffk/UFZygaHjQQMMJilkFfEcEBKEBQuPStWtXVQ44FjE7Es7w6nzEpH9hQU6fPu2zu4kAETjY4DBzPzC/ioUwmAsOR4DM/fffr+bf3Zk8ebJyJqGL6d2NdJ8rDyeIS3Cfp/cFuvvoVmMYNHr0aI+gEsQPYOgxfvx49Rh1g+vBCYjbDfVsgvqGIwtDJnPBDAKO4Cz1tbiooKNy5cpGqVKljCFDhhRY5l9++UWVd8qUKR4LcLAo6L777lP1bJbZHXTXUW7z3jCHXliI5b146PLLL3cNC5OTk9UcP85FjMLf//53j+9ugjiFhx56SH0HnAsnaaSDdCh6P8bZ7mP6SIHxN8af3jcqGhZvqlevbixfvtyw8opFd7CiDQE5BdVjMLMM4SbTxypBgOArCNps2CFIPPbHX4IGA4FC7gE/BYGG6LvvvjO+//57I9LE4X+R60cQcqnLW9CiExJ9KHpCNIPee0I0g6InRDMoekI0wzHeFczLYu4Tc7X+LHghxGmY6bYQdlxYui3HiB6CR7AFIbpz+PDhfGnPHCl6M10TvjAWhNihZ4LoM0SzRTIJot1gvQRfL8ipCMNXVKJWx4je7NJD8HYRPeLpUVbe3L/Begm9Xooa3rJWCdEMip4QzaDoCdEMx4zpSWwxN53wTvwYzNgVa88xfuWY3rNeULfhiJqn6EnIIJmIr2QZoew0i/lmxlt41gtEjwYR8/ChrLen6ElIQKDIPINceObNGIpYzR4DVuVR9Pm998ijiPpGhqVge0IUPQnZyuOGxPxwOJKEUvQF1wvSqyODEvLjo96ReDQYOGgiYYHjb/vUM38pQjSDoidEMyh6oi1Lly5Vm2HAf4AMtNjGOhZguyz3PfqwkUgkoeiJliQnJ6vdjLDBJjYHReppc9/Aojh48KDaiy9cYCcj5P838/mjMYok2nnvM3MuSqkSl7ZaIpHzNGdeuBia9z6v6IUj3pQqXszv9+zatUsJHtt2gSFDhqitriIBtrtyz/tv8vbbb6t8+/jO2GgjWmgj+qNnM2Xsxzvk0KkM+eSx6yQ+ntNBkQKCb/zsyqh/7s6/3yqlS/h3S/dw27wDYKMLzH0jyAh71GHzEWzE8dxzz8l///tf+fLLLwttULD7D3agxcYWADvmQsgnTpxQu+X42gkIG2vgcxB0g/XvsPLYP2/mzJlqF99IoU33vkxignydfEp2H0+TT3dEfg9wYh9ycnLUnvXYQgtbeUHoTz75pNq++tVXX1UbV0Lw2BgT255hhxw0Aua+9HPnzlU9Bmw1hjXt4KOPPlI77VSoUEFtS4ZdbLwPxDVgd55rrrlGNRZoNBB4M2bMmMh+YcMhYIcUfB1zpxRfTFm1x6j95HLjlimfGxcvhr4rTShgAwRsB+3PRghWBptEYHML980isDHE+ewLQR3pWTnG2fQM9W+g7w12p6HRo0erHXjMDSywoUWTJk2MNm3aGEOHDnWdhw09sdEFNgxt27ata2ebc+fOqddxjQULFqi/27dv7/o7ED7//HO1a483+G4oH3bk8a7vQDQAtOneg0Ed68qcr1Jkz4lL1r5b00sbNJLwAqtY2s9uts8xfbxELSJvzZo1ypLDyiLiDeCz4dTDnnOw4ibmhp4A0XDeG1zC2i9evFhuuukm2bp1q/LKFzWmN89x7/KfOnVK7a3ovn9hONGmew8qlCougzrUVX+/+t99kpfHfT50JiUlRe0cC9E3btzY9Ty66NjfHrsXY8NJf7n77rvVLrTYrBIbUpqbmmJMj40svQ9MzeE9GCqYYNdc7BIcKcErDIfgb9fmbEaOcfVzn6pu/n9+uLT/eixwcvc+FMxubDg2BS2MjIwMtVnkgw8+qParMw98Lrr0ffr0UY9r1qxpfPzxx35ft2XLlmpzz4ULF/p1Pvb3a926tRoyYKNS7FP4/PPPR7R7r53orTK2p+hjK/olS5ao+8X72Lhxo1GmTBnj4MGD6rx3331X7RKcnp7u13UnTZqk3n/+/Hm/zsd3HTRokHpPjRo1jLFjx/rctTaconfMXnboksFTij3Bi0qMmZp5QTpOWiNpWbky/Z6W0r1Z9Mf2WJkG7zDGcHZerILlnugm161bN+hVX05ZZbd//37VtYdHfsGCBWG9tlkvOBAc5Ku+/dWAVo4877H9q6v3ybTV++S2q2tw3p6EDEJ50YgvX75crIx9TUwYPPnlSia4PPmEhEp6erry0rs7Ba2ItqKnJ5/oiraiB7T2REe0Fj2tPdERrUUPaO2Jbmgvelp7ohvaix7Q2pNAwHLZb7/9Vs6fPy92hKKntScB8Morr0jDhg1lwIABag38+vXrJVSYLitG0NqTolJiIeIO2XZ27NghO3fulMcff1yeeeYZsVu6LIr+V2jtSVFgueubb76pkmKYS2axDLYocJ6ZcMP9WLZsmcrPZ6bLMp8vU6aMRBKK3g1a+zCB5Rw556N/BLCM5NixYypOHevoAbLldOjQodANIps0aeJa/47xPJbk9u7d25Uuyz3HHnoFiI1HHHxhS2vhGzDTZUHs/fr1U9Y+kmgZe+9PTD6sfdcmjMkPigsZIuNrBvVW1PalVBZB8LejIiX8s5KXuaXFWrRokUqLtW7dOrWGftOmTWpBFDbkhOUFU6dOVdlzwSeffKLEWbt2bVf3Hgk0XnzxRbXoBYtd3NNl4SgIM10W1u9j4RWSdyBd1uuvvy6Rgpa+EGu/YjutvZN57LHHVPccSTLvu+8+tWAG2WxghSHs1q1bu6zyXXfd5XofGgZYb2Dms4NzDwk1zcU2S5YsUQkyigLv/+yzz5TwmzZtKi+//LL8+9//lkhCS+8FV+CFgeKlL1ndaC+txecGQIKPtFhFpcQy33f99dfLtGnT5I9//KOy0nZKl0XR+4C59EIEYvWzm50PjKnjc6GsS9eJIOe80mJ9+OGHhZ7/wQcfyE8//SQjRoxQj7EtN7boNoFlHz9+vM90WQWlwMZ7hg0bJh07dlTPMV1WhDLn+EOks+swc05sM+eAQNNibdmyxShbtqyxePFilQX35ptvNh555BGPc5guy8aij3QuPYo+tqLfGGRaLKS1rlOnjlGxYkVj8ODB+dJiMV2WRdNl+cvUz/aqsX3D6uVkxePh3RWH6bJ8w3RZvmG6rCjBsT0JBKbLcgCM0iOBwHRZDoHz9sRpUPQBWHvM29PaE7tD0fsBY/KJk7Cs6BETjb26rQDH9sRJWFL0P/zwg8yYMUMtgrAKtPaRBVNSp7JS5cf0E3I6+1yhq91IaFhuyg7z1w899JCKh65Xr55YBa7AiwypOenywYHVMnv3x3Iw/Zjr+Tpla8gDjW6Xu+vdJBVKXApnJQ619FhSiCwiWOiAJAM5OTk+z8OCBAQjuB9moxGpY2D72i5r/8m2oyFfD9YskuWN1vFrZGfAx5ojm6X54vvl2c2z5FC6Z+/pUPoJ9Txex3nBfkYgx/nz51VwzenTp6PyeYEcwP3fgn4Lf7BURB7mORFpVKNGDbV66YsvvlA/xOeffy6lSpXyOPf555+XsWPH5rvG3r17pVy5chEr46xvjsrsjcfkyiol5Z17G0t8kFFj+IEQOYUIKjtH5GEhCb4H1pYHEpG37tgWuX/9uEs3sNow1jfxcilv3PzrnpEbLmshkST11++BdfboafoDou+aNWumjI7798fva+a8CweoJyTbQH3/+OOP6r4pXtwz80BaWpo0aNDAXhF5WJYIka9du1aSkpJUyCHWGL/zzjuqy++9Dnn48OGux6j0WrVqSdWqVcMWhuuLYbdUkoVbT0ryqSz5/ue8oKP0zJsC5bV7GC5uNoTN4vC3S//Q15OKFDzA6/GGqPO39J4X1q5+WlqaWglnCrNKlSoyZMgQyczM9PguEBsMkq9kGOZKOtx7JvheSMCxYsUKtcw23OB+QVm9G1l/G11LiR7LFtu1a6cED1DxaEUR0+wNlh76Wn6IComkiCqWSXSN7V9bkyzdmtYMemyPmy3S5Y00KLt7Jld/WJiyRjJzs8UoQvDuwsf5i1LWyoONPNefh0LLli3ll19+8dlYvfbaa66/0QhUrlxZDh8+nO888zsXdJ1wgYbE/KyC7ht/7yNL3W3IE4YKdufQoUOuRIRWgZ780G5eOO3ET8G789buZWH16u/fv1/Onj2rstEiSw7SZEG8OFauXKli6fEcep9YVOQULCV6JBZEamE482D1kZkEGUgwvrcSnLcPHkzHwUsfaI2hV4D3nclJC2t5Dhw4oLrg5lDLpH79+iqZBXoDyKHnbsWRRAMJNJA/Dw0Duv3onZrZbHGt7du3i1WxlOgxTkFusvnz5yuHBObpFy5c6DFesgq09sFxPjcrpPenX/DsCYbqQ/rDH/6gRP/II4+oMb67rwLO49GjR8u9994rGzdudL0GB9o999yjegm+DmS1DccyYy1ED5CGGCmD4AhJTk62TFSeN7T2wVEmITQxlC3uOYsTCg0aNFCeeuSyRzDY5MmTXa9NmjRJ9TpHjhyp7sO2bdsGNHa2ci4ASzny7LzeHivwujfjevuiqJxYXuqUvUzNy/vryANxEie1y9aQSiXCNx179dVXqwO88MILakcbWPzjx4/LnDlzlHVHDjz4mrxjRN5///0C183DB4CeglWxnKW3E1yBFziwgA80Cq73NrjRpT3fIkGbNm1Ugso777xT+vTpI7NmzVI9AW8wB46kl+iJwuGHzSzxLxoKOJ3xN5yNSKZhVSj6EOHYPnAQWlsqIVEF3vgDzsP5fevdGLEyGYahgsIQDYoYioL8SPA3YS4f52MTy3nz5qnnR40aJU899ZQr5z2y26JHYEUo+hDh2D5wEGAzp9PfLs03FyF8MyJvbqenIhKDf/jwYZk5c6aKB0EILmaLsGsN0mJ37txZ/vnPf6rNLmDZERKO+Xs499Dt79Spk3odIHjsjTfeUFN7yGUPZx7SYVsRij4M0NoHTuearWRB5+elZEKiGq/jP3fM5/D6e53Hyg01W4b18zMzM9WuNAi7RT57hHRjXztEgA4dOlQln4TwMX3cokULtcMNpvewJuT+++9X18A+dtiZBhb9qquuUrvYwnsPRx8Ej/3yrIilYu+tlg03kplzmQ33t5DchQfWqMAb71V2gxv1UkOB8sFunFEE33zzjRKx+642voB4IZOaNYPbny8cMBuuBWHm3OBAlx2htYMb9lSBN2k5GVIqrrgkla4U8cbw2muv9eu8ohoFu2FfE2MxdB/bh9phxLgd03m/K1tdKiWWt/Q8dywJR8ecog8jOmbONZd3wtFFIo9Zz97LagOB3ftIZddZvVeLHW/hxUa8OfwToHTp0iFZaTvvcBNJ4APC8l7saIv6dt84M1Ao+giN7feeSNcmSg/z28AUfiiYWWHMJbvEM4kG1qeY9R0sFH2Y0dHaQ5xwdmEmwteWzIEAwcOa4ea286xGuEG9YDoQgg+1MaToI4CO1h6gyxlKt9O8uTFexXQURf8bZu8nHLBWI+3JX71XO08+sTYUfYQ9+aa1J8QqUPQRgivwiFWh6CMIY/KJFaHoI4juUXrEmlD0EYbWnlgNij7C0NoTq0HRRwFae2IlKPooQGtPrARFHyV0XIFHrAlFHyU4b0+sAkUfRTi2J1aAoo+VtV+zX/KckZ6Q2AyKPoYx+ev2n411cYiGUPQxtPazNxzl2J5EHYo+htY++VQWx/Yk6lD0MbL2A9vX+W1sT2tPoghFHyMGdagjZUsUU2N7WnsSTSj6GFG+VHG5u8WlnU0ZpUeiCUUfQ/q1qMZ5exJ1KPoYAsGbY3taexItKHoLjO1p7Uk0oegtMLbnCjwSTSh6C8AVeCSaUPQWgHnySTSh6C0C8+STaEHRWwRaexItKHoLQWtPogFFbyFo7Uk0oOgtBq09iTQUvcWgtSeRhqK3uLVnlB4JNxS9BWGefBJJKHqLwsy5JFJQ9BaF1p5ECorewtDak0hA0VsYWnsSCSh6i8MVeCTcUPQWh3vgkXBD0dsAju1JOKHobQDH9iScUPQ2gdaehAuK3ibQ2pNwQdHbCFp7Eg4oehtBa0/CAUVvMzhvTxwt+q5du8q8efNiXQxLwfX2xLGiX7BggaxcuTLWxbAkzK5DHCf606dPy4gRI6Rhw4axLoolobUnoZAgFgSC7927t2RmZhZ4TnZ2tjpMzp07p/7Ny8tTh9VBGQ3DCLqsA9vXljlfpShr/8m2o9Kt6WXiBEKtF6eS50e9+FtnlhP92rVrZfXq1bJjxw4ZNmxYgedNmDBBxo4dm+/5kydPSlZWllgd/ECpqanqh4yPD67D1feaqjJ74zGZsmq3tKwWL/FxcWJ3wlEvTiTPj3pJS0uzn+gh1ocfflhmzpwp5cqVK/TcMWPGyPDhwz0sfa1ataRq1apSvnx5scOPGBcXp8ob7M097JZKsnDrSTlwKku+/znPEdY+HPXiRPL8qJeSJUvaT/Tjxo2TNm3aSPfu3Ys8NzExUR3eoELscrPgRwylvBXLJKqx/aur98lra5KlW9OaEh9vf2sfar04lbgi6sXf+rKU6N977z3VPa9YsaJ6nJGRIQsXLpRvv/1WZsyYEeviWdaTj7G9GaXnBGtPIoulRL9+/XrJzc11PR45cqS0a9dOBgwYENNy2cGTD2uPKL2uTWo4wtqTyGGp/tMVV1whderUcR1ly5aVpKQkdZCCYUw+sa3ovUE0Hq180TAmnzhG9MR/aO2Jv1D0DoHWnvgLRe8gaO2JP1D0DoLWnvgDRe8waO1JUVD0DoPWnhQFRe9AaO1JYVD0DoTWnhQGRe9QaO1JQVD0DoXWnhQERe9gaO2JLyh6B0NrT3xB0TscWnviDUXvcGjtiTcUvQbQ2hN3KHoNoLUn7lD0msA98IgJRa8J3BWHmFD0GsE98Aig6DWC1p4Ail4zaO0JRa8ZtPaEotcQWnu9oeg1hNZebyh6TXG39ozS0wuKXlMYpacvFL3GMCZfTyh6jaG11xOKXnNo7fWDotccWnv9oOgJV+BpBkVPPKz9tNW09k6HoicKju31gaInCo7t9YGiJy5o7fWAoicuaO31gKInHtDaOx+KnnhAa+98KHqSD87bOxuKnuSD6+2dDUVPfMLsOs6Foic+obV3LhQ9KRBae2dC0ZMCobV3JhQ9KRRae+dB0ZNC4Qo850HRkyJhlJ6zoOhJkTBKz1lQ9MQvaO2dA0VP/ILW3jlQ9MRvaO2dAUVP/IbW3hlQ9CQguALP/lD0JCA4b29/KHoSMBzb2xuKngQMx/b2hqInQUFrb18oehIUtPb2haInQUNrb08sJ/qlS5dKvXr1JCEhQZo3by67du2KdZFIAdDa2xNLiT45OVkGDhwoEydOlCNHjkiDBg1k8ODBsS4WKQRaew1Ev3fv3siURERZdQi+b9++Ur16dRkyZIhs2bIlfB9wMTd81yIKWnv7kRDoG9Dlbtiwodx9993qqFv30g8eDnr06OHxeM+ePVK/fn2f52ZnZ6vD5Ny5c+rfvLw8dfgi7s3rRcpWF6NxL5GG3UVKV5ZYgTIahlFgWe3EwPa1Zc5XKcraf7LtqHRrelnQ13JSvYQTf+rF3zoLWPQnT56UFStWyJIlS2TSpElKlP369ZM+ffpIrVq1JFzk5OTI5MmTZfjw4T5fnzBhgowdO9Zn+bKysvI9Xyz1kFQ9sV3kxHaJS14txn+GS07NdpJ1ZVfJqttFjJKVJJrgB0pNTVU/ZHy8pUZZQdH3mqoye+Mxmbpqt7SsFi/xcXFBXcdp9RIu/KmXtLQ0v64VZ+AqQZKbmytz586Vv/71r8rStm/fXomxY8eOEipjxoxRjcumTZukePHifll6NDpnzpyR8uXL+77oqf0iO5dI3M6lEocG4FeMuGIida+/1ANo1CMqPQD8iGigqlat6oib+1zmBbnu5XWSlpUr/+zfPGhr77R6CRf+1As0UKlSJdU4FKiBYEW/b98++fDDD2Xx4sWyY8cOue2221RXPyMjQ0aPHi3Hjh2TUFizZo3ccccdsmHDBmncuLFf78EXrlChQpFf2MUvaAA+EtmxVOTEtt+eRwNQ73qRJr0j2gDgR/z555+lWrVqjrm5p362V15dvU8aVi8nKx6/TuLjA7f2TqyXcOBPvfirgYC791dffbUcOHBAbr31VnniiSfk9ttvlzJlyqjXUlJSJCkpSUIB1+jfv79Mnz7db8EHRdLvRTqNunScShbZ8ZHqBcjxbSLJay4dH/8lKg2Akzz55tgeK/C6Nwt+bE8iR8CW/t1335VevXpJuXLlwl6YzMxMad26tXTo0EGmTJnieh6NSlwRY8SALX1BmA3AjiWePYD4BDUEkCZ3hKUBcKpFM619g+pl5dPHOwVs7Z1aL1ay9CGN6SMRmINuvS/rX6dOneiI3u8GoFNIPQCn3typmRek46Q1amw//Z6WAVt7p9ZLqDhW9KEQEdEH4gNofIfIVT39bgCcfHOHYu2dXC9WET1rNVAfwJAvRYZ+J3LjMyLVm4oYF38d/z8m8vLvRd7pLfL92yIZp0VXuCuOtaHog24ARl5qAIZ9f6kBqOHWACwbpnUDwD3wrA1FHypVrrzUADxSSAPwj/oi7/xRqwaA1t66UPSRbgAwBMjLFUle7dkAbHlH4rLOiFOhtbcuFH2kGwCXD+BpjwYg/uPHpNrbHSVuwV2O7QG4W3uuwLMOFH2MnIBGjaYSl5er1gF4+AC+m++YBoAr8KwJRR8jJ6Dx0Bdyst9Kyev8tKcPwGGzAFxvbz0o+hhysWIdketGONoJSGtvPSh6yzoBPX0Adm4AaO2tBUVv2QbA3Qdg7waA1t5aUPS2cwL6aABs4ATkHnjWgaK3EzYOBeYeeNaBondyKLDFhgAc21sDit5pTkD3HoDFfAAc21sDit7xPYCCnICx8QHQ2sceil7bWYDY+ABo7WMPRa8LRTkB3XsAEW4AaO1jC0WvI76cgO49gAj7AGjtYwtFrzvuqwGLWg4cxgaA8/axg6InhSwHfiZiTkCut48dFD0peggQISegd3Yd5Gg9lZUqRzJ/Uf86JGer5WA23Bhh26yv3huDhLgzEDLnvrJ2u1T+XYqUqLFTDqb/1tWvU/YyeaBRT7m73k1SoURZ0Zk8psDOD0UfA8KwMcjylG/lgfXjReIuqA1NDPntdoyTS6mzSyUkypxOf5PONVuJruRR9Pmh6GNMEBuDrD36nfxp7fNyUd2CBd+G8ZB/XJws6Py8tsLPo+jzQ9FbCD82Bkm98kZp/ulfJCs3W/IKEby78EsmJMr//jhfy65+Hje7IHZfDfjBvFsl80KWX4IHOC8zN1sWHlgT8eI7HYqeRD0QCElBZ5evUWiXviDe2r2MXv0QoehJ9Pg1DuD0gP/IwRKlxChiJ2Jv4OQ7mH5MzuSkRayIOkDRk6hzPjcrpPenX8gMW1l0hKInUadMQsmQ3l/2QmiNhu5Q9CTqVE4srwJvzHl4f4kzDKmTkymVprWwfE5AK0PRk6iDOXdE2gXxRhlsQPzWzwloZSh6EhMQWotIO8y/+wPOK5VQUvret8w2OQGtCoNzYoSjg3P8xIzIwy2Y50dE3nudx8oNNVv6iARcnD8QKIBQYDvAiDwfUPT2Ff6gL8arwBvgHntv/plYLFHevuHp/IL3NxLQAQ1AHkWfH4revqTmpKtIOwTeYB7epGKxKnL+8FXSqERrWfVYF4mPjwtrKLBc1dM2DQBF7wOK3v6Y6+kPHf9Jate4QhKM0nLdS2slLStXZvyppXRrellwF1YNwJL8i4HiglsOHAsYe08cCcbtmM67vFSS+rdi6RLhyaVnw41BIglFTyxN2DPnVrHHxiCRhKInliaimXOTrL0xSKSg6InliUqe/CrW2xgkUlD0xPJEPU9+knU2BokEFD2xBTHbFScpthuDRAKKntgCS+yKUyU2G4OEG4qe2AZL7YFXJXobg4Qbip7YBktY+xhtDBJOKHpiKyy/B16SmxOwoECgGDcAFD2xFbbaA6+KWyCQuw8gxpGAFD2xHd574NmCKv46ASPfA6Doie2wlbX3ywno5QOI8BCAoie2xJbWPoRAoLgFd0mpXYtEss9JqFD0xJbY3toHshowL1fikldLhc+fFkk/KaFC0RPb4hhr74cTMK/z05L5++6Xng8Rip7YFkdae19A6NeNkNQuUyQcUPTEMdY+5lF6NoGiJ7bGslF6FoaiJ7bHUjH5NoCiJ7aH1j4wKHriCGjt/YeiJ46A1t7Got++fbu0adNGKlWqJKNGjVK50AlxxAo8i2Ap0WdnZ0vPnj2lVatWsnnzZtm5c6fMmzcv1sUiNrT201bT2hdEgliIFStWqN05pkyZIqVLl5bx48fLo48+KgMHDvTZQOBw393D3AkEh9VBGdXGjTYoq53qZWD72jLnqxRl7T/ZdjT4XXFsWC/+1pmlRL9161Zp166dEjxo1qyZsva+mDBhgowdOzbf8ydPnpSsrCyxOviB0MDhh+S2VuGtl77XVJXZG4/J1FW7pWW1eImPC2APPBvXS1pamv1ED2tdt+6l7pm5zVGxYsXkzJkzaozvzpgxY2T48OEe761Vq5ZUrVrVNnvZ4fuhvBR9eOtl2C2VZOHWk5J8Kku+/znPEdY+z496KVmypP1En5CQIImJifm+SEZGRj7R4zzvcwEqxC4iwo9op/LapV4qlklUY/tXV++T19YkS7emNQPb8dam9eJvfVnqbqtcubLqnnt3WUqUKBGzMhF7wnl7m4geU3XffPON63FKSopy1qExICQQOG9vE9F36tRJjc3nzp2rHsN736VLFzWuJyRQOG9vA9FjTP/WW2/J0KFDJSkpSZYuXSqTJk2KdbGITdFmvb2dRQ9uv/12SU5Olvnz58uuXbukcePGsS4SsTGOzq7jFNGDGjVqSPfu3aVKlSqxLgqxObT2NhE9IeGE1t4Tip44Hlp7Tyh6ogW09r9B0RMt4Aq836DoiTYwSu8SFD3RBkbpXYKiJ1oxiNaeoid6UYHWnqIn+jFIc2tP0RPtqKC5tafoiZYM0ngFHkVPtKSCxvP2FD3RlkGaju0peqItFTQd21P0RGsGaWjtKXqiNRU0tPYUPdGeQZpZe4qeaE8Fzaw9RU+I6GXtKXpCRC9rT9ETopm1p+gJ0czaU/SEaGbtKXpCNLP2FD0hmq3Ao+gJ0SxPPkVPiGZ58il6QjSz9hQ9IZpZe4qeEM2sPUVPiGbWnqInRDNrT9ETEoC1d0KUHkVPiGZRehQ9IZrF5FP0hGhm7Sl6QjSz9hQ9IZpZe4qeEM1W4FH0hGi2Bx5FT4hmY3uKnhDNxvYUPSGaWXuKnhDNrD1FT4hm1p6iJ0Qza0/RE6LZvD1FT4hm6+0pekI0y65D0ROimbWn6AnRzNpT9IRoZu0pekI0s/YUPSGarcCj6AnRLEqPoidEsyg9ip4Qzay9pUS/dOlSqVevniQkJEjz5s1l165dsS4SIY6z9pYRfXJysgwcOFAmTpwoR44ckQYNGsjgwYNjXSxCHGftLSN6WHUIvm/fvlK9enUZMmSIbNmyJdbFIsRx1j5BLEKPHj08Hu/Zs0fq169f4PnZ2dnqMDl37pz6Ny8vTx1WB2U0DMMWZY0mTqqXge1ry5yvUpS1/2TbUenW9LKI1ou/dRZ10d9xxx2ybt26fM+/8MILMnToUPV3Tk6OTJ48WYYPH17gdSZMmCBjx47N9/zJkyclKytLrA5+oNTUVPVDxsdbpsMVc5xWL32vqSqzNx6Tqat2S8tq8RIfFxexeklLS/PrWnEGrhJFTpw4IZmZmfmer1y5spQvX179PWbMGFmxYoVs2rRJihcv7relr1Wrlpw5c8Z1HSuDHxENVNWqVR1xc4cLp9XLucwLct3L6yQtK1f+2b950Nben3qBBipVqqQah8I0EHVLj/F6YaxZs0amT58uGzZsKFDwIDExUR3eoELscrPExcXZqrzRwkn1UrFMohrbv7p6n7y2Jlm6Na0p8fFxEakXf+vLUrWakpIi/fv3V6Jv3LhxrItDiCM9+ZYRPbr8cOb16tVLevfuLenp6eqI8uiDEMd78i0j+lWrVsnOnTtl1qxZUq5cOddx6NChWBeNEEdZe8uIHhYeVt37qFOnTqyLRoijrL1lRE+I0xlkEWtP0ROimbWn6AnRzNpT9IRoZu0pekI0s/YUPSGaWXuKnhDN9sCj6AnRLE8+RU+IZnnyKXpCNLP2FD0hmll7ip4Qzaw9RU+IZtaeoidEM2tP0RNiMWsf6Sg9ip4QzaL0KHpCNIvJp+gJ0czaU/SEaGbtKXpCNLP2FD0hmq3Ao+gJsbC1n7Y6/NaeoidEs7E9RU+IZmN7ip4Qzaw9RU+IHcb2a/ZLXpj2daToCbFJTP66/WfDck2KnhCbWPvZG46GZWxP0RNiE2t/4FSWbP0pdGufEJZSEUIiau0n/bGpVErIlha/qxTy9WjpCbEBXa+uIXUrlwrLtSh6QjSDoidEMyh6QjSDoidEMyh6QjSDoidEMyh6QjSDoidEMyh6QjSDoidEMyh6QjSDoidEMyh6QjSDoidEMxyznt74NX/YuXPnxA7k5eVJWlqalCxZUuLj2faasF6Crxfz3je14HjRo0JArVq1Yl0UQmKuhQoVKhT4epxRVLNgo5bw6NGjUq5cOYmLixOrg1YZDdThw4elfPnysS6OZWC9BF8vkDIEX7NmzUJ7SY6x9PiSV1xxhdgN/IC8ufPDegmuXgqz8CYcNBGiGRQ9IZpB0ceIxMREee6559S/5DdYL5GvF8c48ggh/kFLT4hmUPSEaAZFT4hmUPTEMixdulTq1asnCQkJ0rx5c9m1a1esi2Q5unbtKvPmzQvpGhR9DNi+fbu0adNGKlWqJKNGjSoyVloHkpOTZeDAgTJx4kQ5cuSINGjQQAYPHhzrYlmKBQsWyMqVK0O+DkUfZbKzs6Vnz57SqlUr2bx5s+zcuTPkltsJwKpD8H379pXq1avLkCFDZMuWLbEulmU4ffq0jBgxQho2bBjytRwThmsXVqxYIampqTJlyhQpXbq0jB8/Xh599FFl5XSmR48eHo/37Nkj9evXj1l5rAYE37t3b8nMzAz5WrT0UWbr1q3Srl07JXjQrFkzZe3Jb+Tk5MjkyZPlkUceiXVRLMHatWtl9erV8tJLL4XlehR9DFZL1a1b1/UYKwKLFSsmZ86ciWm5rAQiz8qUKcMxvYhkZWXJww8/LDNnzlQrSMMBu/dRBp5p71BKJEbIyMhQjj3dWbNmjUyfPl02bNggxYsXF90ZN26ccvp27949bNek6KNM5cqVlffeHayBLlGihOhOSkqK9O/fX4m+cePGsS6OJXjvvffk5MmTUrFiRfUYxmHhwoXy7bffyowZM4K6JkUfZdBqz5o1y+NGh0cfjYHOwEEFZ16vXr2Uwyo9PV09j26+HZKiRIr169dLbm6u6/HIkSOVT2jAgAFBX5Nj+ijTqVMnNa6fO3euegzvfZcuXdS4XmdWrVqlHJpoEDF2NY9Dhw6JzlxxxRVSp04d11G2bFlJSkpSR7BwlV0MWLZsmerGlipVSmX8WbduHbuzJGpQ9DHi+PHj8t1336muWpUqVWJdHKIRFD0hmsExPSGaQdETohkUPSGaQdETohkUPSGaQdETohkUPSGaQdETohkUPQmZY8eOqY0TsRzWXA/foUMH5v6zKIzII2EB6b+QzXbRokUqqSXWEyCjLbEeFD0JC1j+CZEjDVjbtm3ltddei3WRSAGwe0/ClhEI6a02bdrE3HYWh5aehAXkCMDy4CZNmqg13x9++GGsi0QKgJaehIWnnnpK2rdvr8QOh97y5ctjXSRSALT0JGSQr+3GG2+UHTt2SO3atdVOLGgE8Bjproi1oOgJ0Qx27wnRDIqeEM2g6AnRDIqeEM2g6AnRDIqeEM2g6AnRDIqeEM2g6AnRDIqeENGL/wcBAbIdLV3HOAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "# 1. 生成 x 轴的数据点（从 -1 到 4，生成 50 个均匀分布的点）\n",
    "t = np.linspace(-1, 4, 50)\n",
    "\n",
    "# 2. 创建画布，设置尺寸为 6x5\n",
    "plt.figure(figsize=(6, 5))\n",
    "\n",
    "# 3. 绘制第一条直线：2x + y = 5，转化为 y = 5 - 2x\n",
    "plt.plot(t, 5 - 2*t, label='2x+y=5')\n",
    "\n",
    "# 4. 绘制第二条直线：x + 3y = 5，转化为 y = (5 - x) / 3\n",
    "plt.plot(t, (5 - t)/3, label='x+3y=5')\n",
    "\n",
    "# 5. 绘制交点（即前面求得的解 x[0]=2, x[1]=1）\n",
    "# 颜色设为 '#16a34a' (绿色)，标记大小为 10，标签为 '交点'\n",
    "plt.plot(x[0], x[1], 'o', color='#16a34a', markersize=10, label='交点')\n",
    "\n",
    "# 6. 添加图例、网格、坐标轴标签和标题\n",
    "plt.legend()                     # 显示图例\n",
    "plt.grid(True, alpha=0.3)        # 显示网格，透明度0.3\n",
    "plt.xlabel('x')                  # x轴标签\n",
    "plt.ylabel('y')                  # y轴标签\n",
    "plt.title('唯一解：两条直线相交')   # 标题\n",
    "\n",
    "# 7. 设置坐标轴比例一致，保证图形不会变形\n",
    "plt.gca().set_aspect('equal')\n",
    "\n",
    "# 8. 显示图像\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0576bd2b",
   "metadata": {},
   "source": [
    "手写高斯消元法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "12452891",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "手写消元 x = [ 2.  3. -1.]\n",
      "numpy.solve x = [ 2.  3. -1.]\n",
      "上三角增广矩阵\n",
      "[[ -3.      -1.       2.     -11.    ]\n",
      " [  0.       1.6667   0.6667   4.3333]\n",
      " [  0.       0.       0.2     -0.2   ]]\n"
     ]
    }
   ],
   "source": [
    "def gauss_solve(A, b):\n",
    "    n = len(b)\n",
    "    # 构建增广矩阵 M = [A | b]\n",
    "    M = np.hstack([A.astype(float).copy(), b.reshape(-1, 1).astype(float)])\n",
    "\n",
    "    # 1. 前向消元 (Forward Elimination)\n",
    "    for k in range(n):\n",
    "        # 列主元：寻找当前列下方绝对值最大的元素所在行，防止除以极小数\n",
    "        i = k + int(np.argmax(np.abs(M[k:, k])))\n",
    "        M[[k, i]] = M[[i, k]]  # 交换行\n",
    "\n",
    "        # 如果主元极其微小(接近于0)，说明矩阵奇异，无唯一解\n",
    "        if abs(M[k, k]) < 1e-12:\n",
    "            return None, M\n",
    "\n",
    "        # 执行消元，将主元下方元素化为0\n",
    "        for i in range(k + 1, n):\n",
    "            f = M[i, k] / M[k, k]\n",
    "            M[i] -= f * M[k]\n",
    "\n",
    "    # 2. 回代求解 (Back Substitution)\n",
    "    x = np.zeros(n)\n",
    "    for i in range(n - 1, -1, -1):\n",
    "        # 利用向量化操作 @ 计算已知部分的和\n",
    "        x[i] = (M[i, -1] - M[i, i+1:n] @ x[i+1:n]) / M[i, i]\n",
    "\n",
    "    return x, M\n",
    "\n",
    "# --- 测试第一章的代码 ---\n",
    "A = np.array([[2.0, 1.0, -1.0], [-3.0, -1.0, 2.0], [-2.0, 1.0, 2.0]])\n",
    "b = np.array([8.0, -11.0, -3.0])\n",
    "\n",
    "x, U = gauss_solve(A, b)\n",
    "print('手写消元 x =', x)\n",
    "print('numpy.solve x =', np.linalg.solve(A, b))\n",
    "print('上三角增广矩阵')\n",
    "print(np.round(U, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b8404da",
   "metadata": {},
   "source": [
    "解的情况"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e00984f0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "==== 看下相交 (唯一解) ====\n",
      "增广矩阵消元后:\n",
      "[[2.  1.  5. ]\n",
      " [0.  2.5 2.5]]\n",
      "唯一解 [2. 1.]\n",
      "\n",
      "==== 平行无解 ====\n",
      "增广矩阵消元后:\n",
      "[[ 2.   2.   3. ]\n",
      " [ 0.   0.  -0.5]]\n",
      "无解\n",
      "\n",
      "==== 重合无穷 (无穷多解) ====\n",
      "增广矩阵消元后:\n",
      "[[2. 2. 2.]\n",
      " [0. 0. 0.]]\n",
      "无穷多解或需继续看秩\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# (前提：需先运行第一章的 gauss_solve 函数)\n",
    "# from your_module import gauss_solve\n",
    "\n",
    "def diagnose(A, b):\n",
    "    x, M = gauss_solve(A, b)\n",
    "    print('增广矩阵消元后:')\n",
    "    print(np.round(M, 4))\n",
    "\n",
    "    if x is None:\n",
    "        row = M[-1]\n",
    "        # 判断是否出现 0 = 非零 的矛盾情况 (如 0x + 0y = 3)\n",
    "        if abs(row[-1]) > 1e-8:\n",
    "            print('无解')\n",
    "        else:\n",
    "            print('无穷多解或需继续看秩')\n",
    "    else:\n",
    "        print('唯一解', x)\n",
    "\n",
    "# --- 测试：三种情况的诊断 ---\n",
    "print('==== 看下相交 (唯一解) ====')\n",
    "diagnose(np.array([[2., 1.], [1., 3.]]), np.array([5., 5.]))\n",
    "\n",
    "print('\\n==== 平行无解 ====')\n",
    "diagnose(np.array([[1., 1.], [2., 2.]]), np.array([1., 3.]))\n",
    "\n",
    "print('\\n==== 重合无穷 (无穷多解) ====')\n",
    "diagnose(np.array([[1., 1.], [2., 2.]]), np.array([1., 2.]))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc94513b",
   "metadata": {},
   "source": [
    "病态：扰动0.05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "359df441",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "========== 病态矩阵的扰动分析 ==========\n",
      "原解  [1. 1.]\n",
      "扰动后  [-4.  6.]\n",
      "解的改变量  7.071067811865507\n",
      "cond(A) = 402.0075124842952\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "print(\"========== 病态矩阵的扰动分析 ==========\")\n",
    "\n",
    "# 1. 定义一个病态矩阵 (两行非常接近，行列式接近于0)\n",
    "A = np.array([[1.0, 1.0],\n",
    "              [1.0, 1.01]])\n",
    "\n",
    "# 2. 原始常数向量 b0 与 扰动后的 b1 (仅改变 0.05)\n",
    "b0 = np.array([2.0, 2.01])\n",
    "b1 = np.array([2.0, 2.06])\n",
    "\n",
    "# 3. 分别求解\n",
    "x0 = np.linalg.solve(A, b0)\n",
    "x1 = np.linalg.solve(A, b1)\n",
    "\n",
    "# 4. 打印对比结果与条件数\n",
    "print('原解 ', x0)\n",
    "print('扰动后 ', x1)\n",
    "print('解的改变量 ', np.linalg.norm(x1 - x0))\n",
    "print('cond(A) =', np.linalg.cond(A))\n",
    "\n",
    "# 解释:\n",
    "# 输入仅改变了 0.05，但解的变化量达到了 7.07。\n",
    "# 原因在于 cond(A) 高达 402，说明输入误差被放大了 400 多倍。\n",
    "# 在几何上，这意味着两条直线的斜率极其接近，微小的平移就会让交点位置发生剧烈跳动。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "996eed07",
   "metadata": {},
   "source": [
    "牛顿步就是解方程组"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "775d28c1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "梯度 g = [ 3.2 24. ]\n",
      "牛顿步 d = [-1.6 -1.2]\n",
      "下一步 xy+d = [0. 0.] （二次型应该到原点）\n",
      "错误示范：不要 np.linalg.inv(H) @ (-g)，能 solve 就 solve。\n"
     ]
    }
   ],
   "source": [
    "xy = np.array([1.6, 1.2])\n",
    "g = np.array([2.0 * xy[0], 20.0 * xy[1]])\n",
    "H = np.array([[2.0, 0.0], [0.0, 20.0]])\n",
    "d = np.linalg.solve(H, -g)\n",
    "print('梯度 g =', g)\n",
    "print('牛顿步 d =', d)\n",
    "print('下一步 xy+d =', xy + d, '（二次型应该到原点）')\n",
    "print('错误示范：不要 np.linalg.inv(H) @ (-g)，能 solve 就 solve。')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}