Update eigen slides and exercise

AUTHORS

@@ -1 +1,2 @@
 Rupert Nash
+Joseph Lee

exercises/eigen/explicit.cpp

@@ -6,40 +6,43 @@
 
 int main()
 {
-  int n = 12;
+    int n = 20;
+    int steps = 200;
+    std::vector<double> Avec(n * n);
+
+
+    // Set up matrix A
+    Eigen::Map<Eigen::MatrixXd> A(Avec.data(), n, n);
+    A = Eigen::MatrixXd::Identity(n, n);
+    double delta = 0.4;
+    for (int i = 0; i < n - 1; ++i)
+    {
+        A(i + 1, i) += delta;
+        A(i + 1, i + 1) += -delta;
+
+        A(i, i) += -delta;
+        A(i, i + 1) += +delta;
+    }
+
+    std::cout << "A = \n" << A << std::endl
+              << std::endl;
+
+
+    // T_n
+    Eigen::VectorXd b(n);
+    b.setZero();
+    b.head(n / 2).array() = 1.0;
+
+    std::ofstream f;
+    f.open("explicit_sim.txt");
+    for (int i = 0; i < steps; ++i)
+    {
+        f << b.transpose() << std::endl;
+        // update time-step T_{n+1} = A.T_{n}
+        b = A * b;
+    }
 
-  std::vector<double> Avec(n * n);
-
-  Eigen::Map<Eigen::MatrixXd> A(Avec.data(), n, n);
-  A = Eigen::MatrixXd::Identity(n, n);
-
-  double delta = 0.4;
-
-  for (int i = 0; i < n - 1; ++i)
-  {
-    A(i + 1, i) += delta;
-    A(i + 1, i + 1) += -delta;
-
-    A(i, i) += -delta;
-    A(i, i + 1) += +delta;
-  }
-
-  std::cout << A << std::endl << std::endl;
-
-  Eigen::VectorXd b(n);
-  b.setZero();
-  b.head(n / 2).array() = 1.0;
-
-  std::ofstream f;
-  f.open("a.txt");
-  for (int i = 0; i < 20; ++i)
-  {
     std::cout << b.transpose() << std::endl;
-    f << b << std::endl << std::endl << std::endl;
-    b = A * b;
-  }
-
-  std::cout << b.transpose() << std::endl;
 
-  return 0;
-}
+    return 0;
+}

exercises/eigen/implicit.cpp

@@ -5,38 +5,42 @@
 
 int main()
 {
-  int n = 12;
+    int n = 20;
+    int steps = 100;
 
-  Eigen::MatrixXd A(n, n);
-  A = Eigen::MatrixXd::Identity(n, n);
+    // Set up matrix A
+    Eigen::MatrixXd A(n, n);
+    A = Eigen::MatrixXd::Identity(n, n);
 
-  double delta = -0.4;
+    double delta = -0.4;
 
-  for (int i = 0; i < n - 1; ++i)
-  {
-    A(i + 1, i) += delta;
-    A(i + 1, i + 1) += -delta;
+    for (int i = 0; i < n - 1; ++i)
+    {
+        A(i + 1, i) += delta;
+        A(i + 1, i + 1) += -delta;
 
-    A(i, i) += -delta;
-    A(i, i + 1) += +delta;
-  }
+        A(i, i) += -delta;
+        A(i, i + 1) += +delta;
+    }
 
-  std::cout << A << std::endl << std::endl;
+    std::cout << "A = \n"
+              << A << std::endl
+              << std::endl;
 
-  Eigen::VectorXd b(n);
-  b.setZero();
+    Eigen::VectorXd b(n);
+    b.setZero();
 
-  b.head(n / 2).array() = 1.0;
+    b.head(n / 2).array() = 1.0;
 
-  std::ofstream f;
-  f.open("a.txt");
-  for (int i = 0; i < 20; ++i)
-  {
-    std::cout << b.transpose() << std::endl;
-    f << b << std::endl << std::endl << std::endl;
-    b = A.colPivHouseholderQr().solve(b);
-  }
+    std::ofstream f;
+    f.open("implicit_sim.txt");
+    for (int i = 0; i < 200; ++i)
+    {
+        // Solve A.T_{n+1} = T_{n} for T_{n+1}
+        f << b.transpose() << std::endl;
+        b = A.colPivHouseholderQr().solve(b);
+    }
 
-  std::cout << b.transpose() << std::endl;
-  return 0;
-}
+    std::cout << b.transpose() << std::endl;
+    return 0;
+}

exercises/eigen/movie.py

@@ -0,0 +1,28 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.animation import FuncAnimation
+
+n = 20
+steps = 200
+
+fig, ax = plt.subplots()
+xdata = np.arange(0, n)
+ydata = []
+ln, = ax.plot([], [], 'ro')
+
+data = np.loadtxt("implicit_sim.txt")
+def init():
+    ax.set_xlim(0,n)
+    ax.set_ylim(0,1.1)
+    return ln,
+
+def update(frame):
+    ydata = data[frame]
+    ln.set_data(xdata, ydata)
+    return ln,
+
+ani = FuncAnimation(fig, update, frames=np.arange(0, 200),
+                    init_func=init, blit=True)
+
+
+plt.show()

exercises/eigen/sparse.cpp

@@ -1,42 +1,50 @@
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 
+#include <fstream>
 #include <iostream>
+#include <vector>
 
 int main()
 {
-  int n = 10;
+    int n = 20;
+    int steps = 200;
 
-  Eigen::SparseMatrix<double> A;
-  A.resize(n, n);
+    Eigen::SparseMatrix<double> A;
+    A.resize(n, n);
 
-  double delta = 0.4;
+    double delta = 0.4;
 
-  std::vector<Eigen::Triplet<double>> fill;
-  fill.reserve(n * 4);
+    std::vector<Eigen::Triplet<double>> fill;
+    fill.reserve(n * 4);
 
-  for (int i = 0; i < n - 1; ++i)
-  {
-    fill.push_back(Eigen::Triplet<double>(i + 1, i, delta));
-    fill.push_back(Eigen::Triplet<double>(i + 1, i + 1, -delta));
-    fill.push_back(Eigen::Triplet<double>(i, i, 1.0 - delta));
-    fill.push_back(Eigen::Triplet<double>(i, i + 1, delta));
-  }
-  fill.push_back(Eigen::Triplet<double>(n - 1, n - 1, 1.0));
-  A.setFromTriplets(fill.begin(), fill.end());
+    for (int i = 0; i < n - 1; ++i)
+    {
+        fill.push_back(Eigen::Triplet<double>(i + 1, i, delta));
+        fill.push_back(Eigen::Triplet<double>(i + 1, i + 1, -delta));
+        fill.push_back(Eigen::Triplet<double>(i, i, 1.0 - delta));
+        fill.push_back(Eigen::Triplet<double>(i, i + 1, delta));
+    }
+    fill.push_back(Eigen::Triplet<double>(n - 1, n - 1, 1.0));
+    A.setFromTriplets(fill.begin(), fill.end());
 
-  std::cout << A << std::endl;
+    std::cout << A << std::endl;
 
-  Eigen::VectorXd b(n);
-  b.head(n / 2).array() = 1.0;
-  b.tail(n / 2).array() = 0.0;
+    Eigen::VectorXd b(n);
+    b.head(n / 2).array() = 1.0;
+    b.tail(n / 2).array() = 0.0;
 
-  std::cout << b.transpose() << std::endl;
+    std::cout << b.transpose() << std::endl;
 
-  for (int i = 0; i < 100; ++i)
-    b = A * b;
+    std::ofstream f;
+    f.open("sparse_sim.txt");
+    for (int i = 0; i < steps; ++i)
+    {
+        f << b.transpose() << std::endl;
+        b = A * b;
+    }
 
-  std::cout << b.transpose() << std::endl;
+    std::cout << b.transpose() << std::endl;
 
-  return 0;
-}
+    return 0;
+}