diff --git a/src/simulation/Integrators.ts b/src/simulation/Integrators.ts
index 33324467dd2eec4b1cd375ed47623cb68e7dd748..571b7a4579c674d8e98742ba01f5b96e04493d9d 100644
--- a/src/simulation/Integrators.ts
+++ b/src/simulation/Integrators.ts
@@ -1,3 +1,13 @@
+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());
+}
diff --git a/tests/Integrators.test.ts b/tests/Integrators.test.ts
index 94050ccde9d37c70a73eb2b8cbdefaf29a8fe7d6..cf075b6114b76876c4a14f7e80de917afa387db9 100644
--- a/tests/Integrators.test.ts
+++ b/tests/Integrators.test.ts
@@ -1,41 +1,23 @@
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);
-}
-
-// 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 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);
-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 = "";
for (let i = 0; i < table.length; i++) {
for (let j = 0; j < table[i].length; j++) {
@@ -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", () => {
-
- let k = 0.5
- let c = 2
- let g = (x: number): number => c * math.exp(-k * x);
- let init_value = g(0);
-
- // find the next y value for function g(x)
- let f = (x: Array<number>): Array<number> => [-k * x[0]];
-
- // compare values betwenn [0, 1]
- let bins = 100;
- let h = 1 / bins;
-
- // 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);
- }
+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);
+};
+
+// 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
+
+// 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> = [];
+
+ 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) {
+ let tmp = integrator([result.estimate[result.estimate.length - 1]], h, f);
+ result.estimate.push(tmp[0]);
+ h_list.push(h_accum);
+ h_accum += h;
+ }
+
+ 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);
- // 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]);
+ result.error.forEach(e => {
+ expect(e).toBeLessThan(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);
}
- return result;
+
+ 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++;
+ });
+
+ expect(error_correct_estimate / error_list.length).toBeGreaterThan(error_acc);
+
}
- let max_error = 0.1;
- // 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);
+
+// ------------------------------------------------------------------------------------------------
+
+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("Adaptive step size", () => {
+ adaptive_step_size(Integrator.Midpoint, 0.1, f, "Midpoint-Adaptive");
});
+});
- 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);
+describe("Runge-Kutta", () => {
+ it("Fixed step size", () => {
+ fixed_step_size(Integrator.Runge_Kutta, 0.1, f, "Runge-Kutta-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.Runge_Kutta, 0.1, f, "Runge-Kutta-Adaptive");
});
-});
\ No newline at end of file
+})
\ No newline at end of file
diff --git a/tsconfig.json b/tsconfig.json
new file mode 100644
index 0000000000000000000000000000000000000000..590e673751256a24277dda9877d9f62701ef69ca
--- /dev/null
+++ b/tsconfig.json
@@ -0,0 +1,20 @@
+{
+ "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"]
+}