BiCG_and_BiCGstab.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ax=b\n",
      "行列 A\n",
      "[[1.  2.  0.  0.  0.  0.  0.  0.  0.  0. ]\n",
      " [0.5 1.  2.  0.  0.  0.  0.  0.  0.  0. ]\n",
      " [0.  0.5 1.  2.  0.  0.  0.  0.  0.  0. ]\n",
      " [0.  0.  0.5 1.  2.  0.  0.  0.  0.  0. ]\n",
      " [0.  0.  0.  0.5 1.  2.  0.  0.  0.  0. ]\n",
      " [0.  0.  0.  0.  0.5 1.  2.  0.  0.  0. ]\n",
      " [0.  0.  0.  0.  0.  0.5 1.  2.  0.  0. ]\n",
      " [0.  0.  0.  0.  0.  0.  0.5 1.  2.  0. ]\n",
      " [0.  0.  0.  0.  0.  0.  0.  0.5 1.  2. ]\n",
      " [0.  0.  0.  0.  0.  0.  0.  0.  0.5 1. ]]\n",
      "右辺 b\n",
      "[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
      "解 x\n",
      "[ 2.93000000e+02 -1.46000000e+02  2.50000000e-01  3.68750000e+01\n",
      " -1.80000000e+01  2.81250000e-01  4.85937500e+00 -2.00000000e+00\n",
      "  2.85156250e-01  8.57421875e-01]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "A = np.eye(10) + np.eye(10, k=1)*2 + np.eye(10, k=-1)*0.5\n",
    "b = np.ones(10)\n",
    "x = np.dot(np.linalg.inv(A), b)\n",
    "print(\"Ax=b\")\n",
    "print(\"行列 A\")\n",
    "print(A)\n",
    "print(\"右辺 b\")\n",
    "print(b)\n",
    "print(\"解 x\")\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BiCG Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 -167.21782655596235\n",
      "1000 -5.285692407026592e+17\n",
      "2000 -7.721045015845441e+17\n",
      "3000 nan\n",
      "4000 nan\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:20: RuntimeWarning: divide by zero encountered in double_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5000 nan\n",
      "6000 nan\n",
      "7000 nan\n",
      "8000 nan\n",
      "9000 nan\n",
      "[nan nan nan nan nan nan nan nan nan nan]\n"
     ]
    }
   ],
   "source": [
    "# Initialize solution x\n",
    "x_bicg = np.ones(len(x))\n",
    "A_temp = A.copy()\n",
    "b_base = b.copy()\n",
    "b_temp = np.dot(A_temp, x_bicg)\n",
    "r = b_base - b_temp\n",
    "p = r.copy()\n",
    "ps = r.copy()\n",
    "rs = r.copy()\n",
    "r0 = np.sqrt(np.dot(r, r))\n",
    "eps = 10 ** -10\n",
    "for i in range(10000):\n",
    "    # q _k= Ap_k\n",
    "    b_temp = np.dot(A_temp, p)\n",
    "    # q* _k= AHp*_k\n",
    "    b_temp2 = np.dot(np.linalg.inv(A_temp), ps)\n",
    "    # (r*_k,r_k)\n",
    "    t = np.dot(rs, r)\n",
    "    # alpha_k = (r*_k,r_k) / (p*_k,q_k)\n",
    "    alpha = t / np.dot(ps, b_temp)\n",
    "    # x_k1 = x_k + alpha_k * p_k\n",
    "    x_bicg += alpha * p\n",
    "    # r_k1 = r_k - alpha_k * q_k\n",
    "    r -= alpha * b_temp\n",
    "    # r*_k1 = r*_k - alpha_k * q*_k\n",
    "    rs -= alpha * b_temp2\n",
    "    # (r*_k1, r_k1) / (r*_k, r_k)\n",
    "    beta = np.dot(rs, r) / t\n",
    "    # p_k1 = r_k1 + beta_k * p_k\n",
    "    p = r + beta * p\n",
    "    # p*_k1 = r*_k1 + beta_k * p*_k\n",
    "    ps = rs + beta * ps\n",
    "    if i % 1000 == 0:\n",
    "        print(i, np.sum(x_bicg - x))\n",
    "    if np.sqrt(np.dot(r, r)) / r0 <= eps:\n",
    "        print(i, x_bicg - x)\n",
    "        break\n",
    "print(np.dot(A, x_bicg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [BiCGstab method](https://qiita.com/fukuroder/items/4b708524783192fc2018)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 -167.20930435246737\n",
      "18 [ 8.82954510e-09 -4.41056613e-09 -1.08182074e-11  1.11250387e-09\n",
      " -5.53971091e-10  1.10272902e-12  1.27274191e-10 -6.17257356e-11\n",
      "  3.28246319e-11 -2.61596300e-11]\n",
      "residual:  9.025101595728043e-11\n",
      "[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n"
     ]
    }
   ],
   "source": [
    "# Initialize solution x\n",
    "x_bicgstab = np.ones(len(x))\n",
    "A_temp = A.copy()\n",
    "b_base = b.copy()\n",
    "b_temp = np.dot(A_temp, x_bicgstab)\n",
    "r = b_base - b_temp\n",
    "p = r.copy()\n",
    "rs = r.copy()\n",
    "r0 = np.sqrt(np.dot(r, r))\n",
    "eps = 10 ** -10\n",
    "for i in range(10000):\n",
    "    # q _k= Ap_k\n",
    "    b_temp = np.dot(A_temp, p)\n",
    "    # (r*_k,r_k)\n",
    "    t = np.dot(rs, r)\n",
    "    # alpha_k = (r*_k,r_k) / (p*_k,q_k)\n",
    "    alpha = t / np.dot(rs, b_temp)\n",
    "    # s_k = r_k - alpha_k * q_k\n",
    "    s = r - alpha * b_temp\n",
    "    # w_k = (As_k, s_k) / (As_k, As_k)\n",
    "    As = np.dot(A, s)\n",
    "    w = np.dot(As, s) / np.dot(As, As)\n",
    "    # x_k1 = x_k + alpha_k * p_k\n",
    "    x_bicgstab += alpha * p + w*s\n",
    "    # r_k1 = s_k - w_k * As_k\n",
    "    r = s - w * As\n",
    "    # (r*_k, r_k1) / (r*_k, r_k) * alpha_k / w_k\n",
    "    beta = np.dot(rs, r) / t * alpha / w\n",
    "    # p_k1 = r_k1 + beta_k * (p_k - w_k * Ap_k)\n",
    "    p = r + beta * (p - w * np.dot(A, p)) \n",
    "    if i % 1000 == 0:\n",
    "        print(i, np.sum(x_bicgstab - x))\n",
    "    if np.sqrt(np.dot(r, r)) / r0 <= eps:\n",
    "        print(i, x_bicgstab - x)\n",
    "        print(\"residual: \", np.sqrt(np.dot(r, r)))\n",
    "        break\n",
    "print(np.dot(A, x_bicgstab))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}