Merge pull request #4 from andrei-punko/hyperbolic-eqn-solution

Add test for Hyperbolic equation solution. Compare it with numeric solution

README.MD

@@ -40,4 +40,6 @@ In Russian - algorithm has name "Progonka" (метод прогонки).
 
 ## Usage notes
 
-Check [ParabolicEquationTest](src/test/java/by/andd3dfx/math/pde/equation/ParabolicEquationTest.java)
+Check solved problems in tests:
+- [ParabolicEquationTest](src/test/java/by/andd3dfx/math/pde/equation/ParabolicEquationTest.java)
+- [HyperbolicEquationTest](src/test/java/by/andd3dfx/math/pde/equation/HyperbolicEquationTest.java)

src/test/java/by/andd3dfx/math/pde/equation/HyperbolicEquationTest.java

@@ -0,0 +1,120 @@
+package by.andd3dfx.math.pde.equation;
+
+import by.andd3dfx.math.pde.border.BorderConditionType1;
+import by.andd3dfx.util.FileUtil;
+import org.junit.jupiter.api.Test;
+
+import static java.lang.Math.PI;
+import static java.lang.Math.cos;
+import static java.lang.Math.sin;
+import static org.assertj.core.api.Assertions.assertThat;
+
+/**
+ * Solution of wave equation: Utt = c^2*Uxx
+ * <p>
+ * Check <a href="https://math.libretexts.org/Bookshelves/Differential_Equations/Differential_Equations_(Chasnov)/09%3A_Partial_Differential_Equations/9.06%3A_Solution_of_the_Wave_Equation">article</a>
+ */
+class HyperbolicEquationTest {
+
+    private final double U_MAX = 0.005;     // Max displacement, m
+    private final double L = 0.100;         // Length of string, m
+    private final double TIME = 25;         // Investigated time, sec
+    private final double C_coeff = 1e-2;    // c^2 = T/Ro
+
+    private final double h = L / 2000.0;
+    private final double tau = TIME / 1000.0;
+
+    // We allow difference between numeric & analytic solution no more than 3% of max displacement value
+    private final double EPSILON = U_MAX / 33.;
+
+    @Test
+    void solve() {
+        var waveEquation = buildHyperbolicEquation();
+
+        // Solve equation
+        waveEquation.solve(h, tau);
+
+        // Save numeric solution to file
+        var numericU = waveEquation.gUt(TIME);
+        FileUtil.save(numericU, "./build/hyperbolic-numeric.txt", true);
+
+        // Save analytic solution to file
+        FileUtil.saveFunc(waveEquation.area.x(), (x) -> analyticSolution(x, TIME), "./build/hyperbolic-analytic.txt");
+
+        // Compare numeric & analytic solutions
+        for (var i = 0; i < numericU.getN(); i++) {
+            var x = numericU.x(i);
+            var numericY = numericU.y(i);
+            var analyticY = analyticSolution(x, TIME);
+            assertThat(Math.abs(numericY - analyticY)).isLessThanOrEqualTo(EPSILON);
+        }
+    }
+
+    private HyperbolicEquation buildHyperbolicEquation() {
+        var leftBorderCondition = new BorderConditionType1();
+        var rightBorderCondition = new BorderConditionType1();
+
+        return new HyperbolicEquation(0, L, TIME, leftBorderCondition, rightBorderCondition) {
+            @Override
+            protected double gK(double x, double t, double U) {
+                return C_coeff * C_coeff;
+            }
+
+            @Override
+            protected double gU0(double x) {
+                return getU0(x);
+            }
+        };
+    }
+
+    /**
+     * Initial displacement profile
+     *
+     * @param x space coordinate
+     * @return displacement U_0(x)
+     */
+    private double getU0(double x) {
+        x /= L;
+        if (0 <= x && x <= 0.2) {
+            return U_MAX * 5 * x;
+        }
+        return U_MAX * 1.25 * (1 - x);
+    }
+
+    /**
+     * <pre>
+     * Analytic solution:
+     * U(x,t) = Sum(b_n*u_n(x,t)) = Sum(b_n*sin(n*PI*x/L)*cos(n*PI*C*t/L))
+     *
+     * According to <a href="https://math.libretexts.org/Bookshelves/Differential_Equations/Differential_Equations_(Chasnov)/09%3A_Partial_Differential_Equations/9.06%3A_Solution_of_the_Wave_Equation">article</a>
+     * </pre>
+     *
+     * @param x space coordinate
+     * @param t time
+     * @return concentration C(x,t)
+     */
+    private double analyticSolution(double x, double t) {
+        var result = 0d;
+        for (int n = 1; n <= 100; n++) {
+            result += b(n) * sin(n * PI * x / L) * cos(n * PI * C_coeff * t / L);
+        }
+        return result;
+    }
+
+    /**
+     * Coefficient b_n of a Fourier sine series
+     *
+     * @param n number of coefficient
+     * @return b_n value
+     */
+    private double b(int n) {
+        var integral = 0d;
+        var N = 100d;
+        var dx = L / N;
+        for (int i = 0; i < N; i++) {
+            var x = dx * (i + 0.5);
+            integral += getU0(x) * sin(n * PI * x / L) * dx;
+        }
+        return 2. / L * integral;
+    }
+}

src/test/java/by/andd3dfx/math/pde/equation/ParabolicEquationTest.java

@@ -20,16 +20,16 @@
  */
 class ParabolicEquationTest {
 
-    private final double C0 = 100.0;
-    private final double L = 1e-3;          // 1mm
-    private final double TIME = 1;          // 1sec
-    private final double D = 1e-9;
+    private final double C_MAX = 100.0;
+    private final double L = 0.001;         // Thickness of plate, m
+    private final double TIME = 1;          // Investigated time, sec
+    private final double D = 1e-9;          // Diffusion coefficient
 
     private final double h = L / 100.0;
     private final double tau = TIME / 100.0;
 
     // We allow difference between numeric & analytic solution no more than 1% of max concentration value
-    private final double EPSILON = C0 / 100.;
+    private final double EPSILON = C_MAX / 100.;
 
     @Test
     void solve() {
@@ -40,10 +40,10 @@ void solve() {
 
         // Save numeric solution to file
         var numericU = diffusionEquation.gUt(TIME);
-        FileUtil.save(numericU, "./build/res-numeric.txt", true);
+        FileUtil.save(numericU, "./build/parabolic-numeric.txt", true);
 
         // Save analytic solution to file
-        FileUtil.saveFunc(diffusionEquation.area.x(), (x) -> analyticSolution(x, TIME), "./build/res-analytic.txt");
+        FileUtil.saveFunc(diffusionEquation.area.x(), (x) -> analyticSolution(x, TIME), "./build/parabolic-analytic.txt");
 
         // Compare numeric & analytic solutions
         for (var i = 0; i < numericU.getN(); i++) {
@@ -80,10 +80,10 @@ protected double gU0(double x) {
     private double getU0(double x) {
         x /= L;
         if (0.4 <= x && x <= 0.5) {
-            return C0 * (10 * x - 4);
+            return C_MAX * (10 * x - 4);
         }
         if (0.5 <= x && x <= 0.6) {
-            return C0 * (-10 * x + 6);
+            return C_MAX * (-10 * x + 6);
         }
         return 0;
     }
@@ -102,7 +102,7 @@ private double getU0(double x) {
      */
     private double analyticSolution(double x, double t) {
         var result = 0d;
-        for (int n = 0; n < 100; n++) {
+        for (int n = 1; n <= 100; n++) {
             result += b(n) * sin(n * PI * x / L) * exp(-n * n * PI * PI * D * t / (L * L));
         }
         return result;
@@ -119,7 +119,7 @@ private double b(int n) {
         var N = 100d;
         var dx = L / N;
         for (int i = 0; i < N; i++) {
-            var x = dx * i;
+            var x = dx * (i + 0.5);
             integral += getU0(x) * sin(n * PI * x / L) * dx;
         }
         return 2. / L * integral;