Ch2基础数值方法.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>第二章 基础数值计算方法</center>\n",
    "### <center>许文立</center>\n",
    "### <center>SFU/AHU/CIMERS/国民经济工程实验室(北京)</center>\n",
    "## 2.1 概要\n",
    "\n",
    "* 了解计算经济学中基础的一些数值计算技术\n",
    "* 掌握数值计算在Julia中的应用\n",
    "\n",
    "## 2.2 线性代数回顾\n",
    "定义一个线性方程：n个未知数，n个系数\n",
    "$$a_1 x_1 +a_2 x_2 +……+a_n x_n =b$$\n",
    "所有变量都是一次，且没有非线性函数。多个线性方程组成的系统就是线性方程组：一种紧凑的方式表达线性方程组：\n",
    "$$Ax=b$$\n",
    "其中，A是一个$n \\times n$的系数矩阵，元素为$a_{i,j}$，b为n维向量，x为n为向量。\n",
    "\n",
    "两类解法：\n",
    "* 直接方法\n",
    "* 迭代方法\n",
    "\n",
    "### 2.2.1 直接法\n",
    "#### 2.2.1.1 向前/向后替换法\n",
    "如果系数矩阵是下三角阵，那么，可以向前替换，如果是上三角矩阵，就用向后替换。例如，下列上三角阵：\n",
    "$$A=\\left[\\begin{array}{ccc}a_{11} & a_{12} & a_{13} \\\\0 & a_{22} & a_{23} \\\\0 & 0 & a_{33}\\end{array}\\right]$$\n",
    "\n",
    "那么，我们使用向后替换来解这个方程组：\n",
    "\\begin{array}{r}x_{3}=b_{3} / a_{33} \\\\x_{2}=\\left(b_{2}-a_{23} x_{3}\\right) / a_{22} \\\\x_{1}=\\left(b_{1}-a_{12} x_{2}-a_{13} x_{3}\\right) / a_{11}\\end{array}\n",
    "\n",
    "\n",
    "#### 2.2.1.2 LU因式分解\n",
    "LU分解是最广泛使用的方法。这种算法的原理是，将矩阵A分解为一个下三角矩阵L和一个上三角矩阵U的乘积。\n",
    "\n",
    "LU分解由两个步骤组成：\n",
    "* 将矩阵A分解成L和U，A=LU\n",
    "* 用L和U来解方程组：\n",
    "$$Ax=(LU)x=L(Ux)=b$$\n",
    "\n",
    "定义Ux=y，得到：\n",
    "$$Ly=b$$\n",
    "\n",
    "最后得到解：\n",
    "$$Ux=y$$\n",
    "\n",
    "下面，我们来看一个例子。\n",
    "$$A=\\left[\\begin{array}{lll}1 & 2 & 3 \\\\0 & 3 & 1 \\\\1 & 4 & 2\\end{array}\\right]$$\n",
    "$$b=\\left[\\begin{array}{l}3 \\\\4 \\\\5\\end{array}\\right]$$\n",
    "\n",
    "在Julia中，LU因式分解可以利用基本函数库中的线性代数模块：LinearAlgebra.factorize(A)："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Int64,2}:\n",
       " 1  2  3\n",
       " 0  3  1\n",
       " 1  4  2"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A=[1 2 3;0 3 1;1 4 2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LU{Float64,Array{Float64,2}}\n",
       "L factor:\n",
       "3×3 Array{Float64,2}:\n",
       " 1.0  0.0       0.0\n",
       " 0.0  1.0       0.0\n",
       " 1.0  0.666667  1.0\n",
       "U factor:\n",
       "3×3 Array{Float64,2}:\n",
       " 1.0  2.0   3.0\n",
       " 0.0  3.0   1.0\n",
       " 0.0  0.0  -1.66667"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "f=factorize(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " 1.0  0.0       0.0\n",
       " 0.0  1.0       0.0\n",
       " 1.0  0.666667  1.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L=f.L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " 1.0  2.0   3.0\n",
       " 0.0  3.0   1.0\n",
       " 0.0  0.0  -1.66667"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U=f.U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Int64,1}:\n",
       " 1\n",
       " 2\n",
       " 3"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b=[1,2,3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "只要获得下三角矩阵L和上三角矩阵U，就可以利用向前/向后替换来解方程组："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       "  0.6000000000000003\n",
       "  0.7999999999999999\n",
       " -0.4000000000000001"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\\(A,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在使用LU因式分解时，我们需要主要两个问题：第一个问题是速度，第二个问题是舍入误差。\n",
    "* 关于速度，当我们解一个$Ax=b$类型的方程组时，LU因式分解在大多数情况下是解线性方程组最好的方法。然而，在实践中，尤其是当我们面对重复运算时（例如矩阵A保持不变，而向量b不同），一种更直接的方法则更有效，即计算A的逆矩阵，然后用逆矩阵乘以向量b。\n",
    "* 关于舍入误差问题，只要计算机的精度不足以区分两个接近的数，就会导致最终结果出错。这个问题通常使用转置技术来解决。假设存在舍入误差，任何软件中LU因式分解的实际应用都包括一些转置运算。\n",
    "\n",
    "#### 2.2.1.3 QR分解\n",
    "如果$A^TA$是一个对角矩阵，那么A就称为正交矩阵。对于这个特例，我们可以使用QR分解，即：\n",
    "$$A=QR$$\n",
    "\n",
    "其中，Q为正交矩阵，R为上三角矩阵。\n",
    "\n",
    "我们可以使用下列事实来解初始方程组$Ax=b$：\n",
    "$$Q^TAx=Q^Tb$$  \n",
    "\n",
    "上式等价于：\n",
    "$$Q^TQRx=Q^Tb$$\n",
    "\n",
    "进一步等价于：\n",
    "$$DRx=Q^Tb$$\n",
    "其中，$D=Q^TQ$是一个对角矩阵，因为R是一个上三角矩阵，因此，DR也是一个上三角矩阵。那么，只要我们获得矩阵Q和R，这个方程就很容易使用前后替换法来解。\n",
    "\n",
    "Julia中线性代数模型的基本函数qr(A)可以实现QR分解"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}\n",
       "Q factor:\n",
       "3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:\n",
       " -0.707107   0.301511  -0.639602\n",
       "  0.0       -0.904534  -0.426401\n",
       " -0.707107  -0.301511   0.639602\n",
       "R factor:\n",
       "3×3 Array{Float64,2}:\n",
       " -1.41421  -4.24264  -3.53553\n",
       "  0.0      -3.31662  -0.603023\n",
       "  0.0       0.0      -1.066"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F=qr(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:\n",
       " -0.707107   0.301511  -0.639602\n",
       "  0.0       -0.904534  -0.426401\n",
       " -0.707107  -0.301511   0.639602"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " -1.41421  -4.24264  -3.53553\n",
       "  0.0      -3.31662  -0.603023\n",
       "  0.0       0.0      -1.066"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2.1.4 Cholesky分解\n",
    "对于对称正定矩阵，我们可以使用Cholesky分解。\n",
    "\n",
    "我们可以将Cholesky分解写成：\n",
    "$$A=LL^T$$\n",
    "其中，矩阵L为下三角矩阵。只要实施Cholesky分解，我们就可以像LU分解那样解出初始方程组。Cholesky分解比高斯消元快，因此，在正定矩阵情形下，我们应该优先使用Cholesky分解。\n",
    "\n",
    "在Julia中，我们可以使用函数cholesky()来实施Cholesky分解。\n",
    "\n",
    "#### 2.2.1.5 克莱默法则\n",
    "克莱默法则是一种非常直接的解法，它依赖于使用矩阵A和向量b的元素组成的数学公式。一般不使用它，因为克莱默法则非常慢，只有存在封闭形式解时才可优先使用它。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2.2 迭代法\n",
    "当系数矩阵A为稀疏矩阵时，迭代法非常有效。有两种常用的迭代法：\n",
    "* 高斯-雅克比迭代\n",
    "* 高斯-赛德尔方法\n",
    "\n",
    "为了使用迭代法，将原始线性方程组写成：\n",
    "$$Qx=b+(Q-A)X$$\n",
    "Q是可逆矩阵。进步变换上述方程：\n",
    "$$x=Q^(-1)b+(I-Q^(-1)A)x$$\n",
    "\n",
    "利用上式，我们构建一个关于x的差分方程\n",
    "$$x_{k+1}=Q^(-1)b+(I-Q^(-1)A)x_k$$\n",
    "\n",
    "两种迭代方法有一些共同之处：\n",
    "* 基于A来设立Q，对于高斯-雅克比方法，Q是包含A对角元素的对角矩阵；对于高斯-赛德尔方法，Q从A矩阵的上三角矩阵里得到。\n",
    "* 两种方法都是从解得初值开始，初值要么是向量b，要么是0向量。\n",
    "\n",
    "如果A为对角优势矩阵，即对角元素大于其它元素，那么，两种迭代方法都会收敛。\n",
    "\n",
    "我们使用Julia中的Diagnoal()函数来创建A的对角矩阵，注意：该函数对大小写很敏感。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 1.0000011034417895\n",
       " 1.9999995714315693\n",
       " 2.9999998054451513"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A=[4. 1. -1.;2. 7. 1.;1. -3. 12.];\n",
    "d=Diagonal(A);\n",
    "x=[1.,1.,1.];\n",
    "b=[3.,19.,31.];\n",
    "maxit=1000;\n",
    "dx=[0.,0.,0.];\n",
    "\n",
    "d=zeros(size(A,1));\n",
    "diag=Diagonal(A);\n",
    "\n",
    "for iii=1:size(A,1);\n",
    "    d[iii]=diag[iii,iii]\n",
    "end\n",
    "\n",
    "for i=1:maxit\n",
    "    dx=(b-A*x)./d\n",
    "    x=x+dx\n",
    "    if norm(dx)<0.00001\n",
    "        break\n",
    "    end\n",
    "end\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 1.0000000585033955\n",
       " 2.000000011078294\n",
       " 2.9999999978942906"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# using LinearAlgebra #导入线性代数库\n",
    "    Q=tril(A)\n",
    "    for i=1:maxit\n",
    "        dx=\\(Q,b-A*x)\n",
    "        x=x+dx\n",
    "        if norm(dx)<0.00001\n",
    "            break\n",
    "        end\n",
    "    end\n",
    "    x\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "高斯-赛德尔方法包含一个超松弛参数$\\lambda$，迭代时使用的迭代公式为$x=x+\\lambda dx$，$\\lambda \\in (1,2)$，且它可以提高迭代速度。\n",
    "\n",
    "### 2.2.3 特征值和特征向量\n",
    "\n",
    "给定方阵A，非零特征向量v可以使得下式成立：\n",
    "$$Av=\\lambda v$$\n",
    "其中，$\\lambda$为矩阵A的特征值。Julia内置函数eigvals()可以计算矩阵的特征值："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       "  3.7885167480054305\n",
       "  7.858372678229171\n",
       " 11.353110573765393"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# using LinearAlgebra #导入线性代数库\n",
    "\n",
    "eigvals(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "为了提取特征向量，我们可以输入："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " -0.853082  0.0621391  -0.109341\n",
       "  0.44796   0.815362    0.174559\n",
       "  0.267548  0.575607    0.978557"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigvecs(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "此外，使用函数eigen()可以得到特征值和特征向量。这个函数返回一个目标F，其中包含了特征值F.values和特征向量F.vectors："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}\n",
       "values:\n",
       "3-element Array{Float64,1}:\n",
       "  3.7885167480054305\n",
       "  7.858372678229171\n",
       " 11.353110573765393\n",
       "vectors:\n",
       "3×3 Array{Float64,2}:\n",
       " -0.853082  0.0621391  -0.109341\n",
       "  0.44796   0.815362    0.174559\n",
       "  0.267548  0.575607    0.978557"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F=eigen(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       "  3.7885167480054305\n",
       "  7.858372678229171\n",
       " 11.353110573765393"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " -0.853082  0.0621391  -0.109341\n",
       "  0.44796   0.815362    0.174559\n",
       "  0.267548  0.575607    0.978557"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F.vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们还可以使用eigmin和eigmax来获得最小、最大特征值："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.7885167480054305"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigmin(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11.353110573765393"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigmax(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 插值和曲线拟合\n",
    "### 2.3.1 多项式插值\n",
    "给定n+1个不同的点，我们可构造唯一的n阶多项式来穿过所有的n+1个点：\n",
    "* 拉格朗日插值\n",
    "* 牛顿插值\n",
    "\n",
    "#### 2.3.1.1 拉格朗日插值\n",
    "给定两个点$(x_0,y_0)$和$(x_1,y_1)$，我们给出一条曲线来穿过这两个点：\n",
    "$$y=cx+d$$\n",
    "\n",
    "其中，c表示直线的斜率，d表示截距。这两个参数可以表示为：\n",
    "$$c=\\frac{y_1 -y_0}{x_1 -x_0}$$\n",
    "$$d=y_0 -c x_0$$\n",
    "\n",
    "将参数带入上述直线函数，我们可以写出第一个方程：\n",
    "$$y=\\frac{y_1 -y_0}{x_1 -x_0} x+\\frac{x_1 y_0 -x_0 y_1}{x_1 -x_0}$$\n",
    "\n",
    "而拉格朗日插值法的思想是，将上式转化成包含两个给定点的式子：\n",
    "$$P_1(x)=a_0 (x-x_1)+a_1 (x-x_0)$$\n",
    "将给定的两个点$(x_0,y_0)$和$(x_1,y_1)$，代入上式，可以得到：\n",
    "$$P_1(x)=y_0 \\frac{x-x_1}{x_0-x_1}+y_1 \\frac{x-x_0}{x_1-x_0}$$\n",
    "\n",
    "为了对任意点集合和n阶多项式使用拉格朗日插值法，我们可以根据拉格朗日公式得到：\n",
    "$$P_n(x)=\\sum_{i=0}^{n}y_i l_i(x)$$\n",
    "其中，$l_i(x)=\\frac{x -x_0}{x_1 -x_0}\\frac{x -x_1}{x_2 -x_1}...\\frac{x -x_i}{x_i -x_n}$\n",
    "\n",
    "#### 2.3.1.2 牛顿插值\n",
    "\n",
    "\n",
    "### 2.3.2 样条曲线插值\n",
    "样条曲线的优点在于其平滑特性。这是因为样条曲线是经过给定数据点相同次数的分段连续型多项式连接而成的曲线。它的平滑性不仅仅是因为各点连续，而且其导数也连续。例如，三次样条曲线可以写成：\n",
    "$$s_k(x)=a_k+b_k x+c_k x^2 +d_k x^3$$\n",
    "其中，系数$a_k,b_k,c_k,d_k$是表征每个多项式k的常数，系数在每一个多项式间都不相同。\n",
    "\n",
    "n个样条曲线被用来插值n+1个数据点（$x_i,y_i$）。下面，我们来讲解一下计算分段多项式系数的方法。由于每个插值点$k_i$隐含二阶导数的连续性，我们可以得到：\n",
    "$$f^{'}_{i-1,i}(x_i)=f^{''}_{i,i+1}(x_i)=k_i$$\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "### 2.3.3 B样条曲线插值\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.3",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}