Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Working implementation of Arnoldi algorithm #6

Open
wants to merge 10 commits into
base: master
Choose a base branch
from
Open

Working implementation of Arnoldi algorithm #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
7b5e7d9
Working implementation of Arnoldi algorithm
eramirem Mar 7, 2016
23e3526
Merge remote-tracking branch 'upstream/master'
eramirem Mar 26, 2016
f0bcd08
Fixed 4-stage Runge-Kutta
eramirem Mar 26, 2016
6eea00d
Minor fixes to 4-stage Runge-Kutta Method
eramirem Mar 26, 2016
aaa0b8e
Merge branch 'master' of https://github.com/eramirem/numerical-method…
eramirem Mar 26, 2016
b0aa1b7
Merge remote-tracking branch 'upstream/master'
eramirem Mar 27, 2016
ead65e8
Merge remote-tracking branch 'upstream/master'
eramirem Apr 10, 2016
54d40a9
Fixed some typos
eramirem Apr 12, 2016
a62314c
Merge remote-tracking branch 'upstream/master'
eramirem Apr 12, 2016
9525541
Simplify upwind, minor fixes
eramirem Apr 12, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 27 additions & 14 deletions 06_iterative.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
14 changes: 8 additions & 6 deletions 10_hyperbolic-1.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -223,7 +223,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",
Expand Down Expand Up @@ -367,7 +367,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",
"\n",
"so we can conclude that with the left hand side only that we would be consistent but we have an extra term! "
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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."
]
Expand Down Expand Up @@ -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",
Expand Down