Add solution of diffusion equation with type 2 of border consitions. Fix

 #5

README.MD

@@ -40,4 +40,5 @@ and solve it with help of tridiagonal matrix algorithm (or Thomas algorithm) (or
 
 Check solved problems in tests:
 - [Diffusion problem with type 1 of both condition types](src/test/java/by/andd3dfx/math/pde/solver/ParabolicEquationSolver11Test.java)
+- [Diffusion problem with type 2 of both condition types](src/test/java/by/andd3dfx/math/pde/solver/ParabolicEquationSolver22Test.java)
 - [Wave equation (plucked string) with type 1 of both condition types](src/test/java/by/andd3dfx/math/pde/solver/HyperbolicEquationSolver11Test.java)

src/test/java/by/andd3dfx/math/pde/solver/ParabolicEquationSolver22Test.java

@@ -0,0 +1,132 @@
+package by.andd3dfx.math.pde.solver;
+
+import by.andd3dfx.math.pde.border.BorderConditionType2;
+import by.andd3dfx.math.pde.equation.ParabolicEquation;
+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.exp;
+import static org.assertj.core.api.Assertions.assertThat;
+
+/**
+ * <pre>
+ * Test for ParabolicEquationSolver with type 2 of border conditions on both sides.
+ *
+ * Solution of diffusion equation: Ut = D*Uxx
+ * - sealed pipe ends, so no diffusion through the ends
+ * - initial concentration with triangle profile: most mass concentrated in the center (see method getU0(x))
+ * - constant diffusion coefficient
+ * </pre>
+ *
+ * @see <a href="https://math.libretexts.org/Bookshelves/Differential_Equations/Differential_Equations_(Chasnov)/09%3A_Partial_Differential_Equations/9.05%3A_Solution_of_the_Diffusion_Equation">article</a>
+ */
+class ParabolicEquationSolver22Test {
+
+    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 = C_MAX / 100.;
+
+    @Test
+    void solve() {
+        var diffusionEquation = buildParabolicEquation();
+
+        // Solve equation
+        var solution = new ParabolicEquationSolver()
+                .solve(diffusionEquation, h, tau);
+
+        // Save numeric solution to file
+        var numericU = solution.gUt(TIME);
+        FileUtil.save(numericU, "./build/parabolic22-numeric.txt", true);
+
+        // Save analytic solution to file
+        FileUtil.saveFunc(solution.area().x(), (x) -> analyticSolution(x, TIME), "./build/parabolic22-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 ParabolicEquation buildParabolicEquation() {
+        var leftBorderCondition = new BorderConditionType2();
+        var rightBorderCondition = new BorderConditionType2();
+
+        return new ParabolicEquation(0, L, TIME, leftBorderCondition, rightBorderCondition) {
+            @Override
+            public double gK(double x, double t, double U) {
+                return D;
+            }
+
+            @Override
+            public double gU0(double x) {
+                return getU0(x);
+            }
+        };
+    }
+
+    /**
+     * Initial concentration profile
+     *
+     * @param x space coordinate
+     * @return concentration C_0(x)
+     */
+    private double getU0(double x) {
+        x /= L;
+        if (0.4 <= x && x <= 0.5) {
+            return C_MAX * (10 * x - 4);
+        }
+        if (0.5 <= x && x <= 0.6) {
+            return C_MAX * (-10 * x + 6);
+        }
+        return 0;
+    }
+
+    /**
+     * <pre>
+     * Analytic solution:
+     * U(x,t) = a_0/2 + Sum(a_n*u_n(x,t)) = a_0/2 + Sum(a_n*cos(n*PI*x/L))*exp(-(n*PI/L)^2 * D*t)
+     *
+     * According to <a href="https://math.libretexts.org/Bookshelves/Differential_Equations/Differential_Equations_(Chasnov)/09%3A_Partial_Differential_Equations/9.05%3A_Solution_of_the_Diffusion_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 += a_n(n) * cos(n * PI * x / L) * exp(-Math.pow(n * PI / L, 2.) * D * t);
+        }
+        return a_n(0) / 2. + result;
+    }
+
+    /**
+     * Coefficient a_n of a Fourier sine series
+     *
+     * @param n number of coefficient
+     * @return b_n value
+     */
+    private double a_n(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) * cos(n * PI * x / L) * dx;
+        }
+        return 2. / L * integral;
+    }
+}