6 files changed, 192 insertions, 0 deletions
diff --git a/usth/MATH2.2/project/contour.png b/usth/MATH2.2/project/contour.png
new file mode 100644
index 0000000..dfd4dc1
--- /dev/null
+++ b/usth/MATH2.2/project/contour.png
Binary files differdiff --git a/usth/MATH2.2/project/heatrans.m b/usth/MATH2.2/project/heatrans.m
new file mode 100644
index 0000000..d5d2e04
--- /dev/null
+++ b/usth/MATH2.2/project/heatrans.m
@@ -0,0 +1,15 @@
+function T = heatrans (cp, lambda, rho, Tg, Td, T0, L, dx, dt, N)
+  alpha = lambda / rho / cp;
+  beta = alpha * dt / dx^2;
+  M = round (L / dx);
+  side = repelem (beta, M);
+  A = diag (side, -1) + diag (repelem (1 - 2*beta, M + 1)) + diag (side, 1);
+  A(1, :) = A(end, :) = 0;
+  A(1, 1) = A(end, end) = 1;
+
+  T = repelem (T0, M + 1);
+  [T(1) T(end)] = deal (Tg, Td);
+  for k = 2 : N
+    T(:, k) = A * T(:, k - 1);
+  end
+end
diff --git a/usth/MATH2.2/project/mesh.png b/usth/MATH2.2/project/mesh.png
new file mode 100644
index 0000000..75de2d1
--- /dev/null
+++ b/usth/MATH2.2/project/mesh.png
Binary files differdiff --git a/usth/MATH2.2/project/report.pdf b/usth/MATH2.2/project/report.pdf
new file mode 100644
index 0000000..88b7496
--- /dev/null
+++ b/usth/MATH2.2/project/report.pdf
Binary files differdiff --git a/usth/MATH2.2/project/report.tex b/usth/MATH2.2/project/report.tex
new file mode 100644
index 0000000..e1c36bf
--- /dev/null
+++ b/usth/MATH2.2/project/report.tex
@@ -0,0 +1,155 @@
+\documentclass[a4paper,12pt]{article}
+\usepackage[english,vietnamese]{babel}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{graphicx}
+\usepackage{hyperref}
+\usepackage{lmodern}
+\usepackage{mathtools}
+\usepackage{siunitx}
+
+\newcommand{\exercise}[1]{\noindent\textbf{#1.}}
+\newcommand{\sifrac}[2]{\displaystyle\SI[per-mode=fraction]{#1}{#2}}
+
+\title{Numerical Methods: Heat Transfer}
+\author{Nguyễn Gia Phong--BI9-184}
+\date{\dateenglish\today}
+
+\begin{document}
+\maketitle
+
+Given a bar of length $L = \SI{0.4}{\meter}$ consisting of homogeneous
+and isotropic material with the initial temperature of $T_0 = \SI{0}{\celsius}$.
+Suppose it is perfectly insulated with the exception of the ends with
+the temperature of $T_g = \SI{100}{\celsius}$ and $T_d = \SI{50}{\celsius}$.
+The thermal properties of material will be taken constant.
+\begin{itemize}
+  \item Specific heat capacity:
+    $c_p = \sifrac{900}{\joule\per\kilogram\per\celsius}$
+  \item Thermal conductivity:
+    $\lambda = \sifrac{237}{\watt\per\meter\per\celsius}$
+  \item Density:
+    $\rho = \sifrac{2700}{\kilogram\per\cubic\meter}$
+  \item Thermal diffusivity:
+    $\alpha = \dfrac{\lambda}{\rho c_p}
+            = \sifrac{9.753e-5}{\square\meter\per\second}$
+\end{itemize}
+
+The heat transfer in this bar can be described by the following
+partial differential equation
+\[\frac{\partial T}{\partial t}
+  = \alpha\frac{\partial^2 T}{\partial x^2}\tag{$*$}\]
+where the temperature $T$ is a function of postition $x$ and time $t$.
+
+To solve the problem numerically, we devide space and time into equal intervals
+of norms $\Delta x$ and $\Delta t$ respectively and let $M = L/\Delta x$.
+Consequently, the spartial coordinate is defined as $x_i = (i - 1)\Delta x$
+with $i \in [1\,..\,M + 1]$ and the temporal one is $t_n = n\Delta t$
+with $n \in \mathbb N^*$. With these definitions\footnote{I believe $i = 0$
+and $n = 0$ in the assignment papers are typos since then the domain of $x_i$
+would exceed $L$ and $t_0$ would be negative.}, we denote $T_i^n = T(x_i, t_n)$.
+Using numerical methods, we may start approximating the solutions of ($*$).
+
+\begin{enumerate}
+  \item The left-hand side of ($*$) can be approximated as
+    \[\frac{\partial T}{\partial t} = \frac{T_i^{n+1} - T_i^n}{\Delta t}\]
+  \item Similarly, the right-hand side is expressed in the following form
+    \[\alpha\frac{\partial^2 T}{\partial x^2}
+      = \alpha\frac{T_{i+1}^n - 2T_i^n + T_{i-1}^n}{\Delta x^2}\]
+  \item ($*$) is therefore reformulated as
+    \[\frac{T_i^{n+1} - T_i^n}{\Delta t}
+      = \alpha\frac{T_{i+1}^n - 2T_i^n + T_{i-1}^n}{\Delta x^2}\]
+  \item From the formular above
+    and let $\beta = \alpha\dfrac{\Delta t}{\Delta x^2}$,
+    we get
+    \[T_i^{n+1} = T_i^n + \beta\left(T_{i+1}^n - 2T_i^n + T_{i-1}^n\right)\]
+  \item Boundary conditions:
+    \begin{itemize}
+      \item $\forall n \in \mathbb N^*,\; T_1^n = T(0, t_n) = T_g$
+      \item $\forall n \in \mathbb N^*,\; T_{M+1}^n = T(L, t_n) = T_d$
+      \item $\forall i \in [2\,..\,M],\; T_{i}^1 = T(x_i, 0) = T_0$
+    \end{itemize}
+  \item From (4) and (5), the temperature at point $x_i$ of the bar
+    at time $t_n$ is recursively defined as
+    \[T_i^n = \begin{dcases}
+        T_i^{n-1} + \beta\left(T_{i+1}^{n-1} - 2T_i^{n-1} + T_{i-1}^{n-1}\right)
+        &\text{ if }1 < i \le M \land n > 1\\
+        T_g &\text{ if }i = 1\\
+        T_d &\text{ if }i = M + 1\\
+        T_0 &\text{ otherwise}
+      \end{dcases}\]
+    Since the temperature only depends on the values in the past, values within
+    $(i, n) \in [1\,..\,M+1]\times[1\,..\,N]$ with any $N$ of choice could
+    be computed via dynamic programming:
+    \begin{enumerate}
+      \item Create a 2-dimensional dynamic array $T$ with one-based index
+        and size $(M+1)\times 1$
+      \item Initialize $T$ with $T_1^1 = T_g$, $T_{M+1}^1 = T_d$
+        and $T_i^1 = T_0\ \forall i \in [2\,..\,M]$,
+        where $T_i^n$ is element of row $i$ and column $n$
+      \item For $n = 2$ to $N$
+        \begin{itemize}
+          \item Let $T_1^k = T_g$
+          \item For $i = 2$ to $M$, let $T_i^k = T_i^{k-1}
+            + \beta\left(T_{i+1}^{k-1} - 2T_i^{k-1} + T_{i-1}^{k-1}\right)$
+          \item Let $T_{M+1}^k = T_d$
+        \end{itemize}
+      \item Return $T_i^n$
+    \end{enumerate}
+    Each iteration in (c) can be written in matrix notation as
+    $T^k = AT^{k+1}$, where $T_n$ is column $n$ and $A$ is
+    a matrix of size $(M+1)\times(M+1)$
+    \[A =\begin{bmatrix}
+        1 & 0 & 0 & \cdots & 0 & 0 & 0\\
+        \beta & 1-2\beta & \beta & \cdots & 0 & 0 & 0\\
+        0 & \beta & 1-2\beta & \cdots & 0 & 0 & 0\\
+        \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots\\
+        0 & 0 & 0 & \cdots & 1-2\beta & \beta & 0\\
+        0 & 0 & 0 & \cdots & \beta & 1-2\beta & \beta\\
+        0 & 0 & 0 & \cdots & 0 & 0 & 1
+      \end{bmatrix}\]
+  \item Steps (a) to (c) is then implemented in Octave as
+\begin{verbatim}
+function T = heatrans (cp, lambda, rho, Tg, Td, T0, L,
+                       dx, dt, N)
+  alpha = lambda / rho / cp;
+  beta = alpha * dt / dx^2;
+  M = round (L / dx);
+  side = repelem (beta, M);
+  A = (diag (repelem (1 - 2*beta, M + 1))
+       + diag (side, -1) + diag (side, 1));
+  A(1, :) = A(end, :) = 0;
+  A(1, 1) = A(end, end) = 1;
+
+  T = repelem (T0, M + 1);
+  [T(1) T(end)] = deal (Tg, Td);
+  for k = 2 : N
+    T(:, k) = A * T(:, k - 1);
+  end
+end
+\end{verbatim}
+
+    Choosing $\Delta x = \SI{0.01}{\meter}$, $\Delta t = \SI{0.5}{\second}$
+    and $N = 841$, we define
+\begin{verbatim}
+T = heatrans (900, 237, 2700, 100, 50, 0, 0.4,
+              0.01, 0.5, 841);
+\end{verbatim}
+    then the temperature at point $x_i$ at time $t_n$ is \verb|T(i, n)|.
+    \pagebreak
+
+    To visualize the heat transfer process, we use \verb|mesh| to plot
+    a 3D graph:
+
+    \includegraphics[width=0.841\textwidth]{mesh.png}
+
+    The temperature can be shown more intuitively using \verb|contourf|:
+
+    \includegraphics[width=0.841\textwidth]{contour.png}
+
+    The \verb|script| to reproduce these results along with \verb|heatrans.m|
+    bundled with this report and this document itself are all licensed under a
+    \href{http://creativecommons.org/licenses/by-sa/4.0/}{Creative Commons
+      Attribution-ShareAlike 4.0 International License}.
+\end{enumerate}
+\end{document}
\ No newline at end of file
diff --git a/usth/MATH2.2/project/script b/usth/MATH2.2/project/script
new file mode 100755
index 0000000..99b61f7
--- /dev/null
+++ b/usth/MATH2.2/project/script
@@ -0,0 +1,22 @@
+#!/usr/bin/env octave
+[cp lambda rho] = deal (900, 237, 2700);
+[Tg Td T0] = deal (100, 50, 0);
+[L dx dt N] = deal (0.4, 0.01, 0.5, 841);
+T = heatrans (cp, lambda, rho, Tg, Td, T0, L, dx, dt, N);
+[X Y Z] = deal (0:dx:L, 0:dt:420, T');
+
+figure (1);   # mesh
+mesh (X, Y, Z);
+title ("Temperature of the bar during the first 420 seconds");
+xlabel ("x (m)");
+ylabel ("t (s)");
+zlabel ("T (°C)");
+
+figure (2);   # contour
+contourf (X, Y, Z, 0:7:100, "showtext", "on");
+title ("Temperature of the bar during the first 420 seconds");
+xlabel ("x (m)");
+ylabel ("t (s)");
+
+disp ("Press any key to coninue...");
+kbhit;