From 7b5e7d9161c5be01f83021bf4c0ea4487e5bbd1b Mon Sep 17 00:00:00 2001 From: Elizabeth Ramirez Date: Sun, 6 Mar 2016 19:27:58 -0500 Subject: [PATCH 1/5] Working implementation of Arnoldi algorithm --- 06_iterative.ipynb | 41 +++++++++++++++++++++++++++-------------- 1 file changed, 27 insertions(+), 14 deletions(-) diff --git a/06_iterative.ipynb b/06_iterative.ipynb index 5dbcb0c..7fc865d 100644 --- a/06_iterative.ipynb +++ b/06_iterative.ipynb @@ -1930,21 +1930,34 @@ }, "outputs": [], "source": [ - "U = numpy.ones(m)\n", + "\"\"\"\n", + "k: iterations\n", + "m: size of A\n", + "\"\"\"\n", + "\n", + "k = 10\n", + "m = 15\n", + "x = numpy.linspace(0, 1, m)\n", + "f = numpy.sin(x)\n", + "\n", "A = numpy.eye(m)\n", - "r = numpy.empty((m, m))\n", - "Q = numpy.empty((m, m))\n", - "H = numpy.empty((m, m))\n", - "\n", - "r[:, 0] = f - numpy.dot(A, U)\n", - "Q[:, 0] = r[:, 0] / numpy.linalg.norm(r[:, 0], ord=2)\n", - "for k in xrange(m):\n", - " v = numpy.dot(A, Q[:, k])\n", - " for i in xrange(k):\n", - " H[i, k] = numpy.dot(Q[:, i], v)\n", - " v = v - H[i, k] * Q[:, i] # Orthogonalize\n", - " H[k + 1, k] = numpy.linalg.norm(v, ord=2)\n", - " Q[:, k+1] = v / H[k+1, k]\n", + "H = numpy.zeros((k + 1, k))\n", + "Q = numpy.zeros((m, k + 1))\n", + "U = numpy.ones(m)\n", + "\n", + "r = f - numpy.dot(A, U)\n", + "Q[:, 0] = r / numpy.linalg.norm(r, ord=2)\n", + "\n", + "for j in xrange(k):\n", + " v = numpy.dot(A, Q[:, j])\n", + " for i in xrange(j):\n", + " H[i, j] = numpy.dot(Q[:, i], v)\n", + " v = v - numpy.dot(H[i, j], Q[:, i])\n", + " H[j + 1, j] = numpy.linalg.norm(v, ord=2)\n", + " Q[:, j + 1] = v / H[j + 1, j]\n", + "\n", + "for q in Q.T:\n", + " print numpy.linalg.norm(q, ord=2)\n", " \n", " # GMRES - Check residual R[:, k] to see if it is small enough yet\n", " # If it is then halt and compute U" From f0bcd08cb66dd11aa7484aed0b88b754a5654a6b Mon Sep 17 00:00:00 2001 From: Elizabeth Ramirez Date: Sat, 26 Mar 2016 12:24:50 -0400 Subject: [PATCH 2/5] Fixed 4-stage Runge-Kutta --- 07_ivp.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/07_ivp.ipynb b/07_ivp.ipynb index 835bf12..eea14bd 100644 --- a/07_ivp.ipynb +++ b/07_ivp.ipynb @@ -1295,8 +1295,8 @@ " U_2[n+1] = U_2[n] + 0.5 * delta_t * f(t_n, U_2[n])\n", " U_2[n+1] = U_2[n] + delta_t * f(t_n + 0.5 * delta_t, U_2[n+1])\n", " y_1 = U_4[n]\n", - " y_2 = U_4[n] + 0.5 * delta_t * f(t_n + 0.5 * delta_t, y_1)\n", - " y_3 = U_4[n] + 0.5 * delta_t * f(t_n + 0.5, y_2)\n", + " y_2 = U_4[n] + 0.5 * delta_t * f(t_n, y_1)\n", + " y_3 = U_4[n] + 0.5 * delta_t * f(t_n + 0.5 * delta_t, y_2)\n", " y_4 = U_4[n] + delta_t * f(t_n + 0.5 * delta_t, y_3)\n", " U_4[n+1] = U_4[n] + delta_t / 6.0 * (f(t_n, y_1) + 2.0 * f(t_n + 0.5 * delta_t, y_2) + 2.0 * f(t_n + 0.5 * delta_t, y_3) + f(t_n + delta_t, y_4))\n", " \n", From 6eea00d6a40d35074554af58ca8a4566ad74ea50 Mon Sep 17 00:00:00 2001 From: Elizabeth Ramirez Date: Sat, 26 Mar 2016 12:29:54 -0400 Subject: [PATCH 3/5] Minor fixes to 4-stage Runge-Kutta Method --- 07_ivp.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/07_ivp.ipynb b/07_ivp.ipynb index 835bf12..eea14bd 100644 --- a/07_ivp.ipynb +++ b/07_ivp.ipynb @@ -1295,8 +1295,8 @@ " U_2[n+1] = U_2[n] + 0.5 * delta_t * f(t_n, U_2[n])\n", " U_2[n+1] = U_2[n] + delta_t * f(t_n + 0.5 * delta_t, U_2[n+1])\n", " y_1 = U_4[n]\n", - " y_2 = U_4[n] + 0.5 * delta_t * f(t_n + 0.5 * delta_t, y_1)\n", - " y_3 = U_4[n] + 0.5 * delta_t * f(t_n + 0.5, y_2)\n", + " y_2 = U_4[n] + 0.5 * delta_t * f(t_n, y_1)\n", + " y_3 = U_4[n] + 0.5 * delta_t * f(t_n + 0.5 * delta_t, y_2)\n", " y_4 = U_4[n] + delta_t * f(t_n + 0.5 * delta_t, y_3)\n", " U_4[n+1] = U_4[n] + delta_t / 6.0 * (f(t_n, y_1) + 2.0 * f(t_n + 0.5 * delta_t, y_2) + 2.0 * f(t_n + 0.5 * delta_t, y_3) + f(t_n + delta_t, y_4))\n", " \n", From 54d40a9c3e270f521438576bbd9bcf0b167e6288 Mon Sep 17 00:00:00 2001 From: eramirem Date: Mon, 11 Apr 2016 21:51:44 -0400 Subject: [PATCH 4/5] Fixed some typos --- 10_hyperbolic-1.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/10_hyperbolic-1.ipynb b/10_hyperbolic-1.ipynb index 70adbf6..19d9bfe 100644 --- a/10_hyperbolic-1.ipynb +++ b/10_hyperbolic-1.ipynb @@ -157,7 +157,7 @@ "$$\n", " \\left | \\frac{a \\Delta t}{\\Delta x} \\right | \\leq 1\n", "$$\n", - "which, apart from $a$, shows that $\\Delta t$ and $\\Delta x$ may change at the same order! This leads us to the general conclusion that hyperbolic PDEs are in fact less stiff in general that parabolic PDEs due to the derivatives involved. This is also why hyperbolic PDEs can effectively be solved using explicit time stepping discretizations rather than the implicit ones we found useful for parbolic PDEs. " + "which, apart from $a$, shows that $\\Delta t$ and $\\Delta x$ may change at the same order! This leads us to the general conclusion that hyperbolic PDEs are in fact less stiff in general that parabolic PDEs due to the derivatives involved. This is also why hyperbolic PDEs can effectively be solved using explicit time stepping discretizations rather than the implicit ones we found useful for parabolic PDEs. " ] }, { @@ -222,7 +222,7 @@ "source": [ "This leads to the system $U'(t) = A U(t)$ with\n", "$$\n", - " A = -\\frac{a}{2 \\Delta t} \\begin{bmatrix}\n", + " A = -\\frac{a}{2 \\Delta x} \\begin{bmatrix}\n", " 0 & 1 & & & & -1\\\\\n", " -1 & 0 & 1 \\\\\n", " & -1 & 0 & 1 \\\\\n", @@ -363,7 +363,7 @@ "$$\\begin{aligned}\n", " &\\frac{1}{\\Delta t} \\left [ u + \\Delta t u_t + \\frac{\\Delta t^2}{2} u_{tt} + \\frac{\\Delta t^3}{6} u_{ttt} + \\mathcal{O}(\\Delta t^4) - u \\right ] \\\\\n", " &~~~~~~~~+ \\frac{a}{2 \\Delta x}\\left [ u + \\Delta x u_x + \\frac{\\Delta x^2}{2} u_{xx} + \\frac{\\Delta x^3}{6} u_{xxx} -u + \\Delta x u_x - \\frac{\\Delta x^2}{2} u_{xx} + \\frac{\\Delta x^3}{6} u_{xxx}) + \\mathcal{O}(\\Delta x^4) \\right ] \\\\\n", - " &= u_t + a u_x + \\frac{\\Delta t}{2} u_{tt} + \\frac{\\Delta t^2}{6} u_{ttt} + \\frac{\\Delta x^2}{6} u_{xxx} + \\mathcal{O}(\\Delta t^3) + \\mathcal{O}(\\Delta x^3)\n", + " &= u_t + a u_x + \\frac{\\Delta t}{2} u_{tt} + \\frac{\\Delta t^2}{6} u_{ttt} + a \\frac{\\Delta x^2}{6} u_{xxx} + \\mathcal{O}(\\Delta t^3) + \\mathcal{O}(\\Delta x^3)\n", "\\end{aligned}$$\n", "so we can conclude that with the right hand side only that we would be consistent but we have an extra term! " ] From 9525541dacab7af092d29c0a38a9c77ac7df7ef9 Mon Sep 17 00:00:00 2001 From: eramirem Date: Mon, 11 Apr 2016 23:13:15 -0400 Subject: [PATCH 5/5] Simplify upwind, minor fixes --- 10_hyperbolic-1.ipynb | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/10_hyperbolic-1.ipynb b/10_hyperbolic-1.ipynb index df2cf9e..a683a74 100644 --- a/10_hyperbolic-1.ipynb +++ b/10_hyperbolic-1.ipynb @@ -855,8 +855,10 @@ "t_final = 17.0\n", "while t < t_final:\n", " U_new[0] = U[0] - delta_t / delta_x * (U[0] - U[-1])\n", - " U_new[1:-1] = U[1:-1] - delta_t / delta_x * (U[1:-1] - U[:-2])\n", - " U_new[-1] = U[-1] - delta_t / delta_x * (U[-1] - U[-2])\n", + " #U_new[1:-1] = U[1:-1] - delta_t / delta_x * (U[1:-1] - U[:-2])\n", + " #U_new[-1] = U[-1] - delta_t / delta_x * (U[-1] - U[-2])\n", + " #Periodic values at the boundary will wrap in vectorization operation anyways\n", + " U_new[1:] = U[1:] - delta_t / delta_x * (U[1:] - U[:-1])\n", " U = U_new.copy()\n", " t += delta_t\n", " \n", @@ -1040,7 +1042,7 @@ "$$\n", "and for $a < 0$\n", "$$\n", - " U^{n+1}_j = U^n_j - \\frac{a \\Delta t}{2 \\Delta x} (-3 U^n_j + 4 U^n_{j+1} 0 U^n_{j+2}) + \\frac{a^2 \\Delta t^2}{2 \\Delta x^2} (U^n_j - 2 U^n_{j+1} + U^n_{j+2}).\n", + " U^{n+1}_j = U^n_j - \\frac{a \\Delta t}{2 \\Delta x} (-3 U^n_j + 4 U^n_{j+1} - U^n_{j+2}) + \\frac{a^2 \\Delta t^2}{2 \\Delta x^2} (U^n_j - 2 U^n_{j+1} + U^n_{j+2}).\n", "$$\n", "These methods are stable if $0 \\leq \\nu \\leq 2$ and $-2 \\leq \\nu \\leq 0$ respectively." ] @@ -1278,7 +1280,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", - "version": "2.7.6" + "version": "2.7.11" }, "latex_envs": { "bibliofile": "biblio.bib",