feat: correct unit-test for integrators

src/simulation/Integrators.ts

+import * as math from "mathjs";
+
+export function calculate_error(a: Array<number>, b: Array<number>): Array<number> {
+    let error: Array<number> = [];
+    for (let i = 0; i < a.length; i++) {
+      error.push(Math.abs(a[i] - b[i]));
+    }
+    return error;
+  }
+
 export function Euler(x: Array<number>, h: number, f: (x: Array<number>) => Array<number>): Array<number> {
     let y  = f(x);
     let x_new = x.map((_, i) => x[i] + h * y[i]);
@@ -20,4 +30,13 @@ export function Runge_Kutta(x: Array<number>, h: number, f: (x: Array<number>) =
     return x_new;
 }
 
+export function Adaptive_step_size(integrator: (x: Array<number>, h: number, f: (x: Array<number>) => Array<number>) => Array<number>, 
+    x: Array<number>, h: number, f: (x: Array<number>) => Array<number>, error_max: number, error_list?: Array<number>): number {
+    let x_a = integrator(x, h, f);
+    let x_b = integrator(x, h / 2, f);
+    let error = calculate_error(x_a, x_b);
+    let approximated_error = error.reduce((a, b) => Math.max(a, b));
     
+    if (error_list) error_list.push(approximated_error);
+    return h * math.sqrt(error_max / approximated_error.valueOf());
+}

tests/Integrators.test.ts

 import *  as Integrator from "../src/simulation/Integrators";
 import * as math from "mathjs";
 import * as fs from "fs";
-// select a differential equation with a known solution
 
-// format numbers to a fixed number of decimal places
-function format_number(number: number, decimal_places: number): string {
-    return number.toFixed(decimal_places);
-}
+function write_table(result: Array<number>, exact: Array<number>, error: Array<number>, h_list: Array<number>,
+  title: string, decimal_places: number, error_list?: Array<number>): void {
+    let format_number = (number: number, decimal_places: number): string => Number(number).toFixed(decimal_places);
 
-// find minimum for decimal places
-function decimal_places(n: number): number {
-    let decimal_places = 0;
-    while (n < 1) {
-        n = n * 10;
-        decimal_places++;
-    }
-    return decimal_places+1;
-}
-
-function min_elem(arrs: Array<Array<number>>): number {
-    let min = arrs[0][0];
-    for (let i = 0; i < arrs.length; i++) {
-        for (let j = 0; j < arrs[i].length; j++) {
-            if (arrs[i][j] < min) {
-                min = arrs[i][j];
-            }
-        }
-    }
-    return min;
-}
+  let table = [["h", "exact", "result", "error", "error_estimated"]];
 
-function write_table(result: Array<number>, exact: Array<number>, error: Array<number>, title: string, bins: number, decimal_places: number): void {
-  let table = [["x", "exact", "result", "error"]];
   for (let i = 0; i < result.length; i++) {
-    table.push([format_number(i / bins, decimal_places), format_number(exact[i], decimal_places), format_number(result[i], decimal_places), format_number(error[i], decimal_places)]);
+    table.push([
+      format_number(h_list[i], decimal_places),
+      format_number(exact[i], decimal_places),
+      format_number(result[i], decimal_places),
+      format_number(error[i], decimal_places),
+      (error_list ? format_number(error_list[i], decimal_places) : "n/a")
+    ]);
   }
 
   let table_string = "";
@@ -45,7 +27,7 @@ function write_table(result: Array<number>, exact: Array<number>, error: Array<n
     table_string += "\n";
   }
 
-  fs.writeFileSync(`${title}.txt`, table_string);
+  fs.writeFileSync(`data_${title}.txt`, table_string);
 }
 
 function calculate_error(result: Array<number>, exact: Array<number>): Array<number> {
@@ -56,60 +38,131 @@ function calculate_error(result: Array<number>, exact: Array<number>): Array<num
   return error;
 }
 
-describe("Integrators", () => {
+interface Result {
+  exact: Array<number>;
+  estimate: Array<number>;
+  error: Array<number>;
+}
+
+// exact solution
+function g(x: number): number {
+  return 2 * math.exp(-0.5 * x);
+};
 
-  let k = 0.5
-  let c = 2
-  let g = (x: number): number => c * math.exp(-k * x); 
-  let init_value = g(0);
+// differential equation
+function f(x: Array<number>): Array<number> {
+  return [-0.5 * x[0]];
+}
+
+// INTERVAL CHECK
+const x_max = 2
+const error_max = 0.01;
+const error_acc = 0.8; // 80% of error_estimates are less then max_error
 
-  // find the next y value for function g(x)
-  let f = (x: Array<number>): Array<number> => [-k * x[0]];
+// FIXED STEP SIZE
+function fixed_step_size(integrator: (x: Array<number>, h: number, f: (x: Array<number>) => Array<number>) => Array<number>,
+  h: number, f: (x: Array<number>) => Array<number>, name: string): void {
+    let h_list: Array<number> = [];
 
-  // compare values betwenn [0, 1]
-  let bins = 100;
-  let h = 1 / bins;
+    let result: Result = {
+      exact: [],
+      estimate: [],
+      error: []
+    };
 
-  // calculate the exact solution
-  let exact: Array<number> = [];
-  for(let i = 0; i <= bins; i++) {
-    let x = i / bins;
-    let y = g(x);
-    exact.push(y);
+    h_list.push(0);
+    let h_accum = h;
+    result.estimate.push(g(h_list[0]));
+    while(h_accum < x_max) {
+      let tmp = integrator([result.estimate[result.estimate.length - 1]], h, f);
+      result.estimate.push(tmp[0]);
+      h_list.push(h_accum);
+      h_accum += h;
     }
 
-  // evaluation loop
-  function eval_func(integrator: (x: Array<number>, h: number, f: (x: Array<number>) => Array<number>) => Array<number>, init: number): Array<number> {
-    let result: Array<number> = [init];
-    let x = [init]
-    for(let i = 0; i < bins; i++) {
-      x = integrator(x, h, f);
-      result.push(x[0]);
+    for (let i = 0; i < h_list.length; i++) {
+      result.exact.push(g(h_list[i]));
     }
-    return result;
+
+    result.error = calculate_error(result.estimate, result.exact);
+
+    write_table(result.estimate, result.exact, result.error, h_list, name, 10);
+
+    result.error.forEach(e => {
+      expect(e).toBeLessThan(1);
+    });
+
 }
 
-  let max_error = 0.1;
+function adaptive_step_size(integrator: (x: Array<number>, h: number, f: (x: Array<number>) => Array<number>) => Array<number>,
+  h: number, f: (x: Array<number>) => Array<number>, name: string): void {
+    let h_list: Array<number> = [];
+    let error_list: Array<number> = [];
+
+    let result: Result = {
+      exact: [],
+      estimate: [],
+      error: []
+    };
+
+    h_list.push(0);
+    let h_accum = h;
+    result.estimate.push(g(h_list[0]));
+    while(h_accum < x_max) {
+      h = Integrator.Adaptive_step_size(integrator, [result.estimate[result.estimate.length - 1]], h, f, error_max, error_list);
+      result.estimate.push(integrator([result.estimate[result.estimate.length - 1]], h, f)[0]);
+      h_accum += h;
+      h_list.push(h_accum);
+    }
 
-  // calculate the Euler solution
-  it("Euler", () => {
-    let result = eval_func(Integrator.Euler, init_value);
-    let error = calculate_error(result, exact);
-    error.forEach((e, i) => { expect(e).toBeLessThan(max_error); });
-    write_table(result, exact, error, "Euler", bins, 5);
+    for (let i = 0; i < h_list.length; i++) {
+      result.exact.push(g(h_list[i]));
+    }
+
+    result.error = calculate_error(result.estimate, result.exact);
+
+    write_table(result.estimate, result.exact, result.error, h_list, name, 10, error_list);
+
+    let error_correct_estimate = 0;
+    error_list.forEach(e => {
+      if (e < error_max) error_correct_estimate++;
     });
 
-  it("Midpoint", () => {
-    let result = eval_func(Integrator.Midpoint, init_value);
-    let error = calculate_error(result, exact);
-    error.forEach((e, i) => { expect(e).toBeLessThan(max_error); });
-    write_table(result, exact, error, "Midpoint", bins, 10);
+    expect(error_correct_estimate / error_list.length).toBeGreaterThan(error_acc);
+
+  }
+
+
+
+// ------------------------------------------------------------------------------------------------
+
+describe("Euler", () => {
+
+  it("Fixed step size", () => {
+    fixed_step_size(Integrator.Euler, 0.1, f, "Euler-Fixed");
+  });
+
+  it("Adaptive step size", () => {
+    adaptive_step_size(Integrator.Euler, 0.1, f, "Euler-Adaptive");
+  });
+})
+
+describe("Midpoint", () => {
+  it("Fixed step size", () => {
+    fixed_step_size(Integrator.Midpoint, 0.1, f, "Midpoint-Fixed");
   });
 
-  it("Runge-Kutta", () => {
-    let result = eval_func(Integrator.Runge_Kutta, init_value);
-    let error = calculate_error(result, exact);
-    error.forEach((e, i) => { expect(e).toBeLessThan(max_error); });
-    write_table(result, exact, error, "Runge-Kutta", bins, 20);
+  it("Adaptive step size", () => {
+    adaptive_step_size(Integrator.Midpoint, 0.1, f, "Midpoint-Adaptive");
   });
 });
+
+describe("Runge-Kutta", () => {
+  it("Fixed step size", () => {
+    fixed_step_size(Integrator.Runge_Kutta, 0.1, f, "Runge-Kutta-Fixed");
+  });
+
+  it("Adaptive step size", () => {
+    adaptive_step_size(Integrator.Runge_Kutta, 0.1, f, "Runge-Kutta-Adaptive");
+  });
+})
\ No newline at end of file

tsconfig.json

+{
+  "compilerOptions": {
+    "target": "ESNext",
+    "useDefineForClassFields": true,
+    "module": "ESNext",
+    "lib": ["ESNext", "DOM"],
+    "moduleResolution": "Node",
+    "strict": true,
+    "sourceMap": true,
+    "resolveJsonModule": true,
+    "isolatedModules": true,
+    "esModuleInterop": true,
+    "noEmit": true,
+    "noUnusedLocals": true,
+    "noUnusedParameters": true,
+    "noImplicitReturns": true,
+    "skipLibCheck": true
+  },
+  "include": ["src", "tests/Integrators.test.ts"]
+}