Assignment3.ipynb

{
  "cells":[
    {
      "cell_type":"code",
      "source":[
        "import numpy as np\n",
        "import scipy.sparse as sp"
      ],
      "execution_count":28,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        
      ],
      "execution_count":28,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "#5)\n",
        "offsets = [1, 0, -1]\n",
        "n = 1000\n",
        "v = np.ones( (1,n) )\n",
        "data = np.array([[1, 4, 1]]).T @ v\n",
        "A = sp.dia_matrix( (data, offsets), shape=(n,n) )\n",
        "b = n * np.ones(n)\n",
        "b[0]=5; b[-1] = 5\n",
        "x0 = np.zeros(n)\n",
        "\n",
        "# (CGM(A, b, x0, tol))"
      ],
      "execution_count":29,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"markdown",
      "source":[
        "# Computer problem 2.5.1)"
      ],
      "attachments":{
        
      },
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "def jacobi_method(A, b, x_guess, tol):\n",
        "    n = len(b)\n",
        "    x_correct = np.ones(n)\n",
        "    diag = np.diag(A)\n",
        "    r = A - np.diag(diag)  # reminder matrix without diagonal\n",
        "    back_err = np.max(np.abs(x_guess - x_correct))\n",
        "    steps = 0\n",
        "    while back_err > tol:\n",
        "        steps += 1\n",
        "        x_guess = (b-r.dot(x_guess))\/diag\n",
        "        back_err = np.max(np.abs(x_guess - x_correct))\n",
        "\n",
        "    return back_err, steps\n"
      ],
      "execution_count":30,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "offsets = [1.0, 0.0, -1.0]\n",
        "n = 100\n",
        "v = np.ones( (1,n) )\n",
        "data = np.array([[-1, 3, -1]]).T @ v\n",
        "A = sp.dia_matrix( (data, offsets), shape=(n,n) )\n",
        "b = np.ones(n)\n",
        "b[0], b[-1] = 2.0, 2.0\n",
        "\n",
        "x_guess = np.zeros(n)\n",
        "\n",
        "back_err, step = (jacobi_method(A.toarray(), b, x_guess, 1e-6))\n",
        "print(\"Backward error: \", back_err)\n",
        "print(step, \"steps\")"
      ],
      "execution_count":31,
      "outputs":[
        {
          "name":"stdout",
          "text":[
            "Backward error:  6.867614882333228e-07\n",
            "35 steps\n"
          ],
          "output_type":"stream"
        }
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"markdown",
      "source":[
        "# Computer problem 2.6)\n",
        "**1)**"
      ],
      "attachments":{
        
      },
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "#2.5.1) Conjugate Gradient Method\n",
        "\n",
        "def CGM(A, b, x0, tol):\n",
        "    (n, n) = A.shape\n",
        "    x = []; r = []; d=[]\n",
        "    alpha = []; beta = []\n",
        "    x.append(x0)\n",
        "    r.append(b-A@x[0])\n",
        "    d.append(r[0])\n",
        "    for k in range(n):\n",
        "        if np.linalg.norm(r[k]) < tol:\n",
        "            print(f\"Ferdig etter {k} iterasjoner\")\n",
        "            break\n",
        "        alpha.append( (r[k]).dot(r[k])\/ (d[k].dot(A.dot(r[k]))))\n",
        "        x.append(x[k] + alpha[k]*d[k])\n",
        "        r.append(r[k] - alpha[k]*A.dot(d[k]))\n",
        "        beta.append( (r[k+1].dot(r[k+1])) \/ (r[k].dot(r[k])) )\n",
        "        d.append(r[k+1] + beta[k]*d[k])\n",
        "    return x[-1], k, r[-1]"
      ],
      "execution_count":32,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "# 1a)\n",
        "A1 = np.array([[1, 0], [0, 2]])\n",
        "b1 = np.array([2.0, 4])\n",
        "A2 = np.array([[1, 2], [2, 5]])\n",
        "b2= np.array([1, 1])\n",
        "\n",
        "x0 = [1, 0]\n",
        "tol = 2.0**(-52)\n",
        "\n",
        "x1, k1, r1 = CGM(A1, b1, x0, tol)\n",
        "x2, k2, r2 = CGM(A2, b2, x0, tol)\n",
        "\n",
        "print(x1)\n",
        "print(x2)"
      ],
      "execution_count":33,
      "outputs":[
        {
          "name":"stdout",
          "text":[
            "[2. 2.]\n",
            "[ 3. -1.]\n"
          ],
          "output_type":"stream"
        }
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"markdown",
      "source":[
        "**5)**"
      ],
      "attachments":{
        
      },
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "offsets = [1.0, 0.0, -1.0]\n",
        "n = 10000\n",
        "v = np.ones( (1,n) )\n",
        "data = np.array([[-1, 3, -1]]).T @ v\n",
        "A = sp.dia_matrix( (data, offsets), shape=(n,n) )\n",
        "A = A.toarray()\n",
        "\n",
        "for i in range(n):\n",
        "    if A[i][n-(i+1)] == 0:\n",
        "        A[i][n-(i+1)] = 1\/2\n",
        "\n",
        "b = 1.5 * np.ones(n)\n",
        "b[0], b[n-1] = 2.5, 2.5\n",
        "\n",
        "if (n%2) == 0:\n",
        "    b[int(n\/2)] = 1\n",
        "    b[int(n\/2 -1)] = 1\n",
        "else:\n",
        "    b[int(n\/2)] = 1\n",
        "    b[int(n\/2 + 1)] = 1\n",
        "\n",
        "\n",
        "x0 = np.zeros(n)\n",
        "tol = 2.0**(-52)\n",
        "\n",
        "print(A)\n",
        "\n",
        "x, k, r = CGM(A, b, x0, tol)\n",
        "print(f\"Final residual: {r}\")\n",
        "print(f\"Step length: {k}\")"
      ],
      "execution_count":34,
      "outputs":[
        {
          "name":"stdout",
          "text":[
            "[[ 3.  -1.   0.  ...  0.   0.   0.5]\n",
            " [-1.   3.  -1.  ...  0.   0.5  0. ]\n",
            " [ 0.  -1.   3.  ...  0.5  0.   0. ]\n",
            " ...\n",
            " [ 0.   0.   0.5 ...  3.  -1.   0. ]\n",
            " [ 0.   0.5  0.  ... -1.   3.  -1. ]\n",
            " [ 0.5  0.   0.  ...  0.  -1.   3. ]]\n",
            "Ferdig etter 34 iterasjoner\n",
            "Final residual: [8.48279775e-25 1.97931900e-24 3.77013702e-24 ... 3.77013702e-24\n",
            " 1.97931900e-24 8.48279775e-25]\n",
            "Step length: 34\n"
          ],
          "output_type":"stream"
        }
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"markdown",
      "source":[
        "# Cholesky"
      ],
      "attachments":{
        
      },
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "def Cholesky_transform (A) :\n",
        "    shape_of_A = A.shape\n",
        "    R=np.zeros( shape_of_A )\n",
        "    (n,n)=shape_of_A\n",
        "    for k in range(n) :\n",
        "        if A[k,k]<0: raise ValueError ('A is not positivedefinite')\n",
        "        R[k,k] = np.sqrt(A [k ,k])\n",
        "        u=np.array ([A[k, k + 1:]\/R[k,k]])\n",
        "        R[k, k+1:n]=u\n",
        "        A[k+1:n, k+1:n] = A[k+1:n, k+1:n]-u.T@u\n",
        "    return R"
      ],
      "execution_count":35,
      "outputs":[
        
      ],
      "metadata":{
        
      }
    },
    {
      "cell_type":"code",
      "source":[
        "A = np.array([[4, 2], [2, 3]])\n",
        "print(Cholesky_transform(A))\n",
        "B = np.array([[3, 2, 1], [2, 6, 2], [1, 2, 3]])\n",
        "print(Cholesky_transform(B))"
      ],
      "execution_count":36,
      "outputs":[
        {
          "name":"stdout",
          "text":[
            "[[2.         1.        ]\n",
            " [0.         1.41421356]]\n",
            "[[1.73205081 1.15470054 0.57735027]\n",
            " [0.         2.         0.5       ]\n",
            " [0.         0.         1.        ]]\n"
          ],
          "output_type":"stream"
        }
      ],
      "metadata":{
        
      }
    }
  ],
  "metadata":{
    
  },
  "nbformat":4,
  "nbformat_minor":0
}