diff --git a/_build/images/parabolic-explicit-stencil.png b/_build/images/parabolic-explicit-stencil.png new file mode 100644 index 0000000..f933565 Binary files /dev/null and b/_build/images/parabolic-explicit-stencil.png differ diff --git a/_build/pdes/elliptic.html b/_build/pdes/elliptic.html index b9f4f7f..e7127b5 100644 --- a/_build/pdes/elliptic.html +++ b/_build/pdes/elliptic.html @@ -11,7 +11,7 @@ next_page: url: /quizzes/quiz2-IVPs.html suffix: .ipynb -search: u j equation partial x y points boundary begin end frac where grid point delta using heat example align n laplaces solve figure center values solution method us index formula k transfer right value d square recursion well bmatrix row problems temperature style lets conditions mathbf dimensional column gauss seidel nabla equations img src images png alt width px figcaptionfigure figcaption nine text unknowns matrix elliptic pde left step difference h very vector b mapping set into columns m ghost jacobi shows also finite differences central used direction five stencil consider td above major only need array accurate vdots size +search: u j equation x begin end partial y points boundary n frac where k bmatrix using grid point delta mathbf vdots heat solve values b method example align solution laplaces right figure center us value index formula transfer text xi d square recursion well quad row old problems left temperature style lets conditions dimensional column jacobi nabla equations img src images png alt width px figcaptionfigure difference figcaption nine unknowns very matrix cdots elliptic pde step h vector mapping into columns m ghost iterative gauss seidel shows also finite differences central used direction five stencil consider td above system major comment: "***PROGRAMMATICALLY GENERATED, DO NOT EDIT. SEE ORIGINAL FILES IN /content***" --- @@ -141,18 +141,18 @@

Example: heat transfer in a sq
-
clear all; clc
-
-A = [
-1 0 0 0  0 0 0 0 0;
-0 1 0 0  0 0 0 0 0;
-0 0 1 0  0 0 0 0 0;
-0 0 0 1  0 0 0 0 0;
-0 1 0 1 -4 1 0 1 0;
-0 0 0 0  0 1 0 0 0;
-0 0 0 0  0 0 1 0 0;
-0 0 0 0  0 0 0 1 0;
-0 0 0 0  0 0 0 0 1];
+
clear all; clc
+
+A = [
+1 0 0 0  0 0 0 0 0;
+0 1 0 0  0 0 0 0 0;
+0 0 1 0  0 0 0 0 0;
+0 0 0 1  0 0 0 0 0;
+0 1 0 1 -4 1 0 1 0;
+0 0 0 0  0 1 0 0 0;
+0 0 0 0  0 0 1 0 0;
+0 0 0 0  0 0 0 1 0;
+0 0 0 0  0 0 0 0 1];
 b = [100; 100; 100; 100; 0; 100; 100; 0; 100];
 
 % solve system of linear equations
@@ -435,7 +435,7 @@ 
Example: heat transfer
 
 

     

-
clear; clc; close all
+
clear; clc; close all
 
 h = 0.1;
 x = [0 : h : 1]; n = length(x);
@@ -449,43 +449,43 @@ 
Example: heat transfer
 u_left = 100;
 u_right = 100;
 u_bottom = 100;
-u_top = 0;
+u_top = 0;
 
-for j = 1 : m
-    for i = 1 : n
+for j = 1 : m
+    for i = 1 : n
         % for convenience we calculate all the indices once
         kij = (j-1)*n + i;
         kim1j = (j-1)*n + i - 1;
         kip1j = (j-1)*n + i + 1;
         kijm1 = (j-2)*n + i;
-        kijp1 = j*n + i;
-        
-        if i == 1 
+        kijp1 = j*n + i;
+        
+        if i == 1 
             % this is the left boundary
             A(kij, kij) = 1;
-            b(kij) = u_left;
-        elseif i == n 
+            b(kij) = u_left;
+        elseif i == n 
             % right boundary
             A(kij, kij) = 1;
-            b(kij) = u_right;
-        elseif j == 1 
+            b(kij) = u_right;
+        elseif j == 1 
             % bottom boundary
             A(kij, kij) = 1;
-            b(kij) = u_bottom;
-        elseif j == m 
+            b(kij) = u_bottom;
+        elseif j == m 
             % top boundary
             A(kij, kij) = 1;
-            b(kij) = u_top;
-        else
-            % these are the coefficients for the interior points,
+            b(kij) = u_top;
+        else
+            % these are the coefficients for the interior points,
             % based on the recursion formula
             A(kij, kim1j) = 1;
             A(kij, kip1j) = 1;
             A(kij, kijm1) = 1;
             A(kij, kijp1) = 1;
-            A(kij, kij) = -4;
-        end
-    end
+            A(kij, kij) = -4;
+        end
+    end
 end
 u = A \ b;
 
@@ -565,7 +565,7 @@ 
Neumann (derivative) boundary
 
 

     

-
clear; clc; close all
+
clear; clc; close all
 
 h = 0.1;
 x = [0 : h : 1]; n = length(x);
@@ -578,45 +578,45 @@ 
Neumann (derivative) boundary
 
 u_left = 100;
 u_right = 100;
-u_bottom = 0;
+u_bottom = 0;
 
-for j = 1 : m
-    for i = 1 : n
+for j = 1 : m
+    for i = 1 : n
         % for convenience we calculate all the indices once
         kij = (j-1)*n + i;
         kim1j = (j-1)*n + i - 1;
         kip1j = (j-1)*n + i + 1;
         kijm1 = (j-2)*n + i;
-        kijp1 = j*n + i;
-        
-        if i == 1 
+        kijp1 = j*n + i;
+        
+        if i == 1 
             % this is the left boundary
             A(kij, kij) = 1;
-            b(kij) = u_left;
-        elseif i == n 
+            b(kij) = u_left;
+        elseif i == n 
             % right boundary
             A(kij, kij) = 1;
-            b(kij) = u_right;
-        elseif j == 1 
+            b(kij) = u_right;
+        elseif j == 1 
             % bottom boundary
             A(kij, kij) = 1;
-            b(kij) = u_bottom;
-        elseif j == m 
+            b(kij) = u_bottom;
+        elseif j == m 
             % top boundary, using the ghost node + recursion formula
             A(kij, kim1j) = 1;
             A(kij, kip1j) = 1;
             A(kij, kijm1) = 2;
-            A(kij, kij) = -4;
-        else
-            % these are the coefficients for the interior points,
+            A(kij, kij) = -4;
+        else
+            % these are the coefficients for the interior points,
             % based on the recursion formula
             A(kij, kim1j) = 1;
             A(kij, kip1j) = 1;
             A(kij, kijm1) = 1;
             A(kij, kijp1) = 1;
-            A(kij, kij) = -4;
-        end
-    end
+            A(kij, kij) = -4;
+        end
+    end
 end
 u = A \ b;
 
@@ -689,43 +689,43 @@ 
Iterative solutions for (
 u_left = 100;
 u_right = 100;
 u_bottom = 100;
-u_top = 0;
+u_top = 0;
 
-for j = 1 : m
-    for i = 1 : n
+for j = 1 : m
+    for i = 1 : n
         % for convenience we calculate all the indices once
         kij = (j-1)*n + i;
         kim1j = (j-1)*n + i - 1;
         kip1j = (j-1)*n + i + 1;
         kijm1 = (j-2)*n + i;
-        kijp1 = j*n + i;
-        
-        if i == 1 
+        kijp1 = j*n + i;
+        
+        if i == 1 
             % this is the left boundary
             A(kij, kij) = 1;
-            b(kij) = u_left;
-        elseif i == n 
+            b(kij) = u_left;
+        elseif i == n 
             % right boundary
             A(kij, kij) = 1;
-            b(kij) = u_right;
-        elseif j == 1 
+            b(kij) = u_right;
+        elseif j == 1 
             % bottom boundary
             A(kij, kij) = 1;
-            b(kij) = u_bottom;
-        elseif j == m 
+            b(kij) = u_bottom;
+        elseif j == m 
             % top boundary
             A(kij, kij) = 1;
-            b(kij) = u_top;
-        else
-            % these are the coefficients for the interior points,
+            b(kij) = u_top;
+        else
+            % these are the coefficients for the interior points,
             % based on the recursion formula
             A(kij, kim1j) = 1;
             A(kij, kip1j) = 1;
             A(kij, kijm1) = 1;
             A(kij, kijp1) = 1;
-            A(kij, kij) = -4;
-        end
-    end
+            A(kij, kij) = -4;
+        end
+    end
 end
 u = A \ b;
 
@@ -778,21 +778,21 @@ 
Iterative solutions for (
 
 

     

-
clear all; clc;
+
clear all; clc;
 
 step_sizes = [0.1, 0.05, 0.025, 0.02, 0.0125, 0.01];
 
 n = length(step_sizes);
-nums = zeros(n,1); times = zeros(n,1);
+nums = zeros(n,1); times = zeros(n,1);
 
-for i = 1 : n
-    [times(i), nums(i)] = heat_equation(step_sizes(i));
-end
+for i = 1 : n
+    [times(i), nums(i)] = heat_equation(step_sizes(i));
+end
 
-loglog(nums, times, '-o')
+loglog(nums, times, '-o')
 xlabel('Number of unknowns'); 
 ylabel('Time for direct solution (sec)')
-hold on
+hold on
 
 x = nums(3:end);
 n2 = x.^2 * (times(3) / x(1)^2);
@@ -848,7 +848,30 @@ 
Iterative solutions for (
 
 

 

-
Jacobi method 

+
Jacobi method 
The Jacobi method essentially works by starting with an initial guess to the solution, then using the recursion formula to solve for values at each point, then repeating this until the values converge (i.e., stop changing).

+
An algorithm we can use to solve Laplace's equation:

+
    
+
  1. Set some initial guess for all unknowns: $u_{i,j}^{\text{old}}$
    2. 
+
  2. Set the boundary values
    3. 
+
  3. For each point in the interior, use the recursion formula to solve for new values based on old values at the surrounding points: $u_{i,j} = \left( u_{i+1,j}^{\text{old}} + u_{i-1,j}^{\text{old}} + u_{i,j+1}^{\text{old}} + u_{i,j-1}^{\text{old}} \right)/4$.
    4. 
+
  4. Check for convergence: is $\epsilon$ less than some tolerance, such as $10^{-6}$? Where $\epsilon = \max \left| u_{i,j} - u_{i,j}^{\text{old}} \right|$. If no, then return to step 2 and repeat.
    5. 
+

+
More formally, if we have a system $A \mathbf{x} = \mathbf{b}$, where
+\begin{equation}
+A = \begin{bmatrix}
+a_{11} & a_{12} & \cdots & a_{1n} \\
+a_{21} & a_{22} & \cdots & a_{2n} \\
+\vdots & \vdots & \ddots & \vdots \\
+a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}
+\quad \mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix}
+\quad \mathbf{b} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{bmatrix}
+\end{equation}
+then we can solve iterative for $\mathbf{x}$ using
+\begin{equation}
+x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j \neq i} a_{ij} x_j^{(k)} \right) , \quad i = 1,2,\ldots, n 
+\end{equation}
+where $x_i^{(k)}$ is a value of the solution at iteration $k$ and $x_i^{(k+1)}$ is at the next iteration.

+
 

 

 

@@ -886,43 +909,43 @@ 
Jacobi method% dummy value for residual variable
 epsilon = 1.0; 
 
-num_iter = 0;
-while epsilon > 1e-6
+num_iter = 0;
+while epsilon > 1e-6
     u_old = u;
     
-    epsilon = 0;
-    for j = 1 : m
-        for i = 1 : n
+    epsilon = 0;
+    for j = 1 : m
+        for i = 1 : n
             kij = (j-1)*n + i;
             kim1j = (j-1)*n + i - 1;
             kip1j = (j-1)*n + i + 1;
             kijm1 = (j-2)*n + i;
-            kijp1 = j*n + i;
+            kijp1 = j*n + i;
 
-            if i == 1 
+            if i == 1 
                 % this is the left boundary
-                u(kij) = u_left;
-            elseif i == n 
+                u(kij) = u_left;
+            elseif i == n 
                 % right boundary
-                u(kij) = u_right;
-            elseif j == 1 
+                u(kij) = u_right;
+            elseif j == 1 
                 % bottom boundary
-                u(kij) = u_bottom;
-            elseif j == m 
+                u(kij) = u_bottom;
+            elseif j == m 
                 % top boundary
-                u(kij) = u_top;
-            else
-                % interior points
-                u(kij) = (u_old(kip1j) + u_old(kim1j) + u_old(kijm1) + u_old(kijp1))/4.0;
-            end
-        end
-    end
-    
-    epsilon = max(abs(u - u_old));
-    num_iter = num_iter + 1;
-end
-
-u_square = reshape(u, [n, m]);
+                u(kij) = u_top;
+            else
+                % interior points
+                u(kij) = (u_old(kip1j) + u_old(kim1j) + u_old(kijm1) + u_old(kijp1))/4.0;
+            end
+        end
+    end
+    
+    epsilon = max(abs(u - u_old));
+    num_iter = num_iter + 1;
+end
+
+u_square = reshape(u, [n, m]);
 % the "20" indicates the number of levels for the contour plot
 contourf(x, y, u_square', 20);
 c = colorbar;
@@ -964,11 +987,11 @@ 
Jacobi method
step_sizes = [0.1, 0.05, 0.025, 0.02, 0.0125, 0.01, 0.005];
 n = length(step_sizes);
 
-nums_jac = zeros(n,1); times_jac = zeros(n,1); num_iter_jac = zeros(n,1);
+nums_jac = zeros(n,1); times_jac = zeros(n,1); num_iter_jac = zeros(n,1);
 
-for i = 1 : n
-    [times_jac(i), nums_jac(i), num_iter_jac(i)] = heat_equation_jacobi(step_sizes(i));
-end
+for i = 1 : n
+    [times_jac(i), nums_jac(i), num_iter_jac(i)] = heat_equation_jacobi(step_sizes(i));
+end
 

 
     

@@ -1015,15 +1038,22 @@ 
Jacobi method

 

-
Gauss-Seidel method 
The Gauss-Seidel method essentially works by starting with an initial guess to the solution, then using the recursion formula to solve for values at each point, then repeating this until the values converge (i.e., stop changing).

-
An algorithm we can use to solve Laplace's equation:

-
    
-
  1. Set some initial guess for all unknowns: $u_{i,j}^{\text{old}}$
    2. 
-
  2. Set the boundary values, and set $u_{i,j} = u_{i,j}^{\text{old}}$
    3. 
-
  3. For each point in the interior, use the recursion formula to solve for new values based on surrounding points: $u_{i,j} = \left( u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} \right)/4$.
    4. 
-
  4. Check for convergence: is $\epsilon$ less than some tolerance, such as $10^{-6}$? Where $\epsilon = \max \left| u_{i,j} - u_{i,j}^{\text{old}} \right|$. If no, then return to step 2 and repeat.
    5. 
-

-
One important component of the Gauss-Seidel method is to immediately incorporate updated values in the solution; this contrasts the Jacobi method, which only uses old values to calculate new ones. The Gauss-Seidel generally converges faster.

+
Gauss-Seidel method 
The Gauss-Seidel method is very similar to the Jacobi method, but with one important difference: rather than using all old values to calculate the new values, incorporate updated values as they are available. Because the method incorporates newer information more quickly, it tends to converge faster (meaning, with fewer iterations) than the Jacobi method.

+
Formally, if we have a system $A \mathbf{x} = \mathbf{b}$, where
+\begin{equation}
+A = \begin{bmatrix}
+a_{11} & a_{12} & \cdots & a_{1n} \\
+a_{21} & a_{22} & \cdots & a_{2n} \\
+\vdots & \vdots & \ddots & \vdots \\
+a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}
+\quad \mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix}
+\quad \mathbf{b} = \begin{bmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{bmatrix}
+\end{equation}
+then we can solve iterative for $\mathbf{x}$ using
+\begin{equation}
+x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+1)} - \sum_{j =i+1}^n a_{ij} x_j^{(k)} \right) , \quad i = 1,2,\ldots, n 
+\end{equation}
+where $x_i^{(k)}$ is a value of the solution at iteration $k$ and $x_i^{(k+1)}$ is at the next iteration.

 
 

 

@@ -1062,43 +1092,43 @@ 
Gauss-Seidel method% dummy value for residual variable
 epsilon = 1.0; 
 
-num_iter = 0;
-while epsilon > 1e-6
+num_iter = 0;
+while epsilon > 1e-6
     u_old = u;
     
-    epsilon = 0;
-    for j = 1 : m
-        for i = 1 : n
+    epsilon = 0;
+    for j = 1 : m
+        for i = 1 : n
             kij = (j-1)*n + i;
             kim1j = (j-1)*n + i - 1;
             kip1j = (j-1)*n + i + 1;
             kijm1 = (j-2)*n + i;
-            kijp1 = j*n + i;
+            kijp1 = j*n + i;
 
-            if i == 1 
+            if i == 1 
                 % this is the left boundary
-                u(kij) = u_left;
-            elseif i == n 
+                u(kij) = u_left;
+            elseif i == n 
                 % right boundary
-                u(kij) = u_right;
-            elseif j == 1 
+                u(kij) = u_right;
+            elseif j == 1 
                 % bottom boundary
-                u(kij) = u_bottom;
-            elseif j == m 
+                u(kij) = u_bottom;
+            elseif j == m 
                 % top boundary
-                u(kij) = u_top;
-            else
-                % interior points
-                u(kij) = (u(kip1j) + u(kim1j) + u(kijm1) + u(kijp1))/4.0;
-            end
-        end
-    end
-    
-    epsilon = max(abs(u - u_old));
-    num_iter = num_iter + 1;
-end
-
-u_square = reshape(u, [n, m]);
+                u(kij) = u_top;
+            else
+                % interior points
+                u(kij) = (u(kip1j) + u(kim1j) + u(kijm1) + u(kijp1))/4.0;
+            end
+        end
+    end
+    
+    epsilon = max(abs(u - u_old));
+    num_iter = num_iter + 1;
+end
+
+u_square = reshape(u, [n, m]);
 %% the "20" indicates the number of levels for the contour plot
 contourf(x, y, u_square', 20);
 c = colorbar;
@@ -1140,13 +1170,13 @@ 
Gauss-Seidel method
step_sizes = [0.1, 0.05, 0.025, 0.02, 0.0125, 0.01, 0.005];
 n = length(step_sizes);
 
-nums_gs = zeros(n,1); times_gs = zeros(n,1); num_iter_gs = zeros(n,1);
+nums_gs = zeros(n,1); times_gs = zeros(n,1); num_iter_gs = zeros(n,1);
 
-for i = 1 : n
-    [times_gs(i), nums_gs(i), num_iter_gs(i)] = heat_equation_gaussseidel(step_sizes(i));
-end
+for i = 1 : n
+    [times_gs(i), nums_gs(i), num_iter_gs(i)] = heat_equation_gaussseidel(step_sizes(i));
+end
 
-loglog(nums, times, '-o'); hold on
+loglog(nums, times, '-o'); hold on
 loglog(nums_jac, times_jac, '-^')
 loglog(nums_gs, times_gs, '-x')
 xlabel('Number of unknowns'); 
@@ -1202,7 +1232,7 @@ 
Gauss-Seidel method
     

 
loglog(nums_jac, num_iter_jac, '-^')
-hold on;
+hold on;
 loglog(nums_gs, num_iter_gs, '-x')
 xlabel('Number of unknowns')
 ylabel('Number of iterations required')
diff --git a/content/images/parabolic-explicit-stencil.png b/content/images/parabolic-explicit-stencil.png
new file mode 100644
index 0000000..f933565
Binary files /dev/null and b/content/images/parabolic-explicit-stencil.png differ
diff --git a/content/pdes/elliptic.ipynb b/content/pdes/elliptic.ipynb
index 90af4fb..6bd61c9 100644
--- a/content/pdes/elliptic.ipynb
+++ b/content/pdes/elliptic.ipynb
@@ -711,7 +711,31 @@
    "metadata": {},
    "source": [
     "### Jacobi method\n",
-    "\n"
+    "\n",
+    "The Jacobi method essentially works by starting with an initial guess to the solution, then using the recursion formula to solve for values at each point, then repeating this until the values converge (i.e., stop changing). \n",
+    "\n",
+    "An algorithm we can use to solve Laplace's equation:\n",
+    "\n",
+    "1. Set some initial guess for all unknowns: $u_{i,j}^{\\text{old}}$\n",
+    "2. Set the boundary values\n",
+    "3. For each point in the interior, use the recursion formula to solve for new values based on old values at the surrounding points: $u_{i,j} = \\left( u_{i+1,j}^{\\text{old}} + u_{i-1,j}^{\\text{old}} + u_{i,j+1}^{\\text{old}} + u_{i,j-1}^{\\text{old}} \\right)/4$.\n",
+    "4. Check for convergence: is $\\epsilon$ less than some tolerance, such as $10^{-6}$? Where $\\epsilon = \\max \\left| u_{i,j} - u_{i,j}^{\\text{old}} \\right|$. If no, then return to step 2 and repeat.\n",
+    "\n",
+    "More formally, if we have a system $A \\mathbf{x} = \\mathbf{b}$, where\n",
+    "\\begin{equation}\n",
+    "A = \\begin{bmatrix}\n",
+    "a_{11} & a_{12} & \\cdots & a_{1n} \\\\\n",
+    "a_{21} & a_{22} & \\cdots & a_{2n} \\\\\n",
+    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
+    "a_{n1} & a_{n2} & \\cdots & a_{nn} \\end{bmatrix}\n",
+    "\\quad \\mathbf{x} = \\begin{bmatrix} x_1 \\\\ x_2 \\\\ \\vdots \\\\ x_n \\end{bmatrix}\n",
+    "\\quad \\mathbf{b} = \\begin{bmatrix} b_1 \\\\ b_2 \\\\ \\vdots \\\\ b_n \\end{bmatrix}\n",
+    "\\end{equation}\n",
+    "then we can solve iterative for $\\mathbf{x}$ using\n",
+    "\\begin{equation}\n",
+    "x_i^{(k+1)} = \\frac{1}{a_{ii}} \\left( b_i - \\sum_{j \\neq i} a_{ij} x_j^{(k)} \\right) , \\quad i = 1,2,\\ldots, n \n",
+    "\\end{equation}\n",
+    "where $x_i^{(k)}$ is a value of the solution at iteration $k$ and $x_i^{(k+1)}$ is at the next iteration."
    ]
   },
   {
@@ -845,21 +869,28 @@
    "source": [
     "### Gauss-Seidel method\n",
     "\n",
-    "The Gauss-Seidel method essentially works by starting with an initial guess to the solution, then using the recursion formula to solve for values at each point, then repeating this until the values converge (i.e., stop changing). \n",
-    "\n",
-    "An algorithm we can use to solve Laplace's equation:\n",
-    "\n",
-    "1. Set some initial guess for all unknowns: $u_{i,j}^{\\text{old}}$\n",
-    "2. Set the boundary values, and set $u_{i,j} = u_{i,j}^{\\text{old}}$\n",
-    "3. For each point in the interior, use the recursion formula to solve for new values based on surrounding points: $u_{i,j} = \\left( u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} \\right)/4$.\n",
-    "4. Check for convergence: is $\\epsilon$ less than some tolerance, such as $10^{-6}$? Where $\\epsilon = \\max \\left| u_{i,j} - u_{i,j}^{\\text{old}} \\right|$. If no, then return to step 2 and repeat.\n",
+    "The Gauss-Seidel method is very similar to the Jacobi method, but with one important difference: rather than using all old values to calculate the new values, incorporate updated values as they are available. Because the method incorporates newer information more quickly, it tends to converge faster (meaning, with fewer iterations) than the Jacobi method.\n",
     "\n",
-    "One important component of the Gauss-Seidel method is to immediately incorporate updated values in the solution; this contrasts the Jacobi method, which only uses old values to calculate new ones. The Gauss-Seidel generally converges faster."
+    "Formally, if we have a system $A \\mathbf{x} = \\mathbf{b}$, where\n",
+    "\\begin{equation}\n",
+    "A = \\begin{bmatrix}\n",
+    "a_{11} & a_{12} & \\cdots & a_{1n} \\\\\n",
+    "a_{21} & a_{22} & \\cdots & a_{2n} \\\\\n",
+    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
+    "a_{n1} & a_{n2} & \\cdots & a_{nn} \\end{bmatrix}\n",
+    "\\quad \\mathbf{x} = \\begin{bmatrix} x_1 \\\\ x_2 \\\\ \\vdots \\\\ x_n \\end{bmatrix}\n",
+    "\\quad \\mathbf{b} = \\begin{bmatrix} b_1 \\\\ b_2 \\\\ \\vdots \\\\ b_n \\end{bmatrix}\n",
+    "\\end{equation}\n",
+    "then we can solve iterative for $\\mathbf{x}$ using\n",
+    "\\begin{equation}\n",
+    "x_i^{(k+1)} = \\frac{1}{a_{ii}} \\left( b_i - \\sum_{j=1}^{i-1} a_{ij} x_j^{(k+1)} - \\sum_{j =i+1}^n a_{ij} x_j^{(k)} \\right) , \\quad i = 1,2,\\ldots, n \n",
+    "\\end{equation}\n",
+    "where $x_i^{(k)}$ is a value of the solution at iteration $k$ and $x_i^{(k+1)}$ is at the next iteration."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 56,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [
     {
diff --git a/content/pdes/parabolic.ipynb b/content/pdes/parabolic.ipynb
new file mode 100644
index 0000000..8ef8600
--- /dev/null
+++ b/content/pdes/parabolic.ipynb
@@ -0,0 +1,72 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Parabolic PDEs\n",
+    "\n",
+    "A classic example of a parabolic partial differential equation (PDE) is the one-dimensional unsteady heat equation:\n",
+    "\\begin{equation}\n",
+    "\\frac{\\partial T}{\\partial t} = \\alpha \\frac{\\partial^2 T}{\\partial t^2} \n",
+    "\\end{equation}\n",
+    "where $T(x,t)$ is the temperature varying in space and time, and $\\alpha$ is the thermal diffusivity: $\\alpha = k / (\\rho c_p)$, which is a constant.\n",
+    "\n",
+    "We can solve this using finite differences to represent the spatial derivatives and time derivatives separately.\n",
+    "First, let's rearrange the PDE slightly:\n",
+    "\\begin{equation}\n",
+    "\\frac{\\partial^2 T}{\\partial x^2} = \\frac{1}{\\alpha} \\frac{\\partial T}{\\partial t}\n",
+    "\\end{equation}\n",
+    "\n",
+    "Let's use a *central difference* for the spatial derivative with a spacing of $\\Delta x$, and a *forward difference* for the time derivative with a time-step size of $\\Delta t$. With these choices, we can obtain an approximation to the PDE that applies at time $t^k$ and spatial location $x_i$:\n",
+    "\\begin{equation}\n",
+    "\\frac{T_{i-1}^k - 2 T_i^k + T_{i+1}^k}{\\Delta x^2} = \\frac{1}{\\alpha} \\left( \\frac{T_i^{k+1} - T_i^k}{\\Delta t} \\right)\n",
+    "\\end{equation}\n",
+    "where $T_i^k$ is the temperature at time $t^k$ and spatial location $x_i$. The following figure shows the stencil of points involved in the PDE, for a domain with five points in the $x$-direction.\n",
+    "\n",
+    "
\n",
+    "  
\n",
+    "  \"stencil\n",
+    "  
Figure: Stencil for explicit solution to heat equation
\n",
+    "  
\n",
+    "
\n",
+    "\n",
+    "To solve the heat equation for a one-dimensional domain over $0 \\leq x \\leq L$, we will need both initial conditions at $t = 0$ and boundary conditions at $x=0$ and $x=L$ (for all time). In terms of our nodal values, this means we need $T_i^{k=1}$ for $i = 1 \\ldots n$, where $n$ is the number of points, as well as information about $T_1^k$ and $T_n^k$ for all times $k$.\n",
+    "\n",
+    "We can rearrange the above equation to obtain our recursion formula:\n",
+    "\\begin{equation}\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "L"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Matlab",
+   "language": "matlab",
+   "name": "matlab"
+  },
+  "language_info": {
+   "codemirror_mode": "octave",
+   "file_extension": ".m",
+   "help_links": [
+    {
+     "text": "MetaKernel Magics",
+     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
+    }
+   ],
+   "mimetype": "text/x-octave",
+   "name": "matlab",
+   "version": "0.16.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}