diff --git a/tests/Integrators.test.ts b/tests/Integrators.test.ts index cf075b6114b76876c4a14f7e80de917afa387db9..132557db852b90b9b346fa0f3836fca790e229dc 100644 --- a/tests/Integrators.test.ts +++ b/tests/Integrators.test.ts @@ -4,7 +4,7 @@ import * as fs from "fs"; 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); + let format_number = (number: number, decimal_places: number): string => Number(number).toFixed(decimal_places); let table = [["h", "exact", "result", "error", "error_estimated"]]; @@ -30,6 +30,34 @@ function write_table(result: Array<number>, exact: Array<number>, error: Array<n fs.writeFileSync(`data_${title}.txt`, table_string); } +function write_table_2d(result: Array<Array<number>>, exact: Array<Array<number>>, error: Array<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); + + let table = [["h", "exact", "result", "error", "error_estimated"]]; + + for (let i = 0; i < result.length; i++) { + table.push([ + format_number(h_list[i], decimal_places), + exact[i].map(e => format_number(e, decimal_places)).join(" , "), + result[i].map(e => format_number(e, decimal_places)).join(" , "), + error[i].map(e => format_number(e, decimal_places)).join(" , "), + (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++) { + table_string += table[i][j]; + if (j < table[i].length - 1) table_string += " | "; + } + table_string += "\n"; + } + + fs.writeFileSync(`data_${title}.txt`, table_string); +} + function calculate_error(result: Array<number>, exact: Array<number>): Array<number> { let error: Array<number> = []; for (let i = 0; i < result.length; i++) { @@ -38,131 +66,190 @@ function calculate_error(result: Array<number>, exact: Array<number>): Array<num return error; } -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]]; +function calculate_error_2d(result: Array<Array<number>>, exact: Array<Array<number>>): Array<Array<number>> { + let error: Array<Array<number>> = []; + for (let i = 0; i < result.length; i++) { + error.push([]); + for (let j = 0; j < result[i].length; j++) { + error[i].push(Math.abs(result[i][j] - exact[i][j])); + } + } + return error; } // INTERVAL CHECK const x_max = 2 -const error_max = 0.01; +const error_max = 0.1; 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); +function eval_1d_integrator(integrator: (x: Array<number>, h: number, f: (x: Array<number>) => Array<number>) => Array<number>, + h: number, f: (x: Array<number>) => Array<number>, g: (x: number) => Array<number>, name: string, adaptive: boolean = false): void { - result.error.forEach(e => { - expect(e).toBeLessThan(1); - }); + let h_list: Array<number> = []; + let error_list: Array<number> = []; -} + let exact: Array<number> = []; + let estimate: Array<number> = []; + let error: Array<number> = []; -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); - } + h_list.push(0); + let h_accum = h; + estimate.push(g(h_list[0])[0]); + while (h_accum < x_max) { + if (adaptive) { + h = Integrator.Adaptive_step_size(integrator, [estimate[estimate.length - 1]], h, f, error_max, (adaptive) ? error_list : undefined); + } + estimate.push(integrator([estimate[estimate.length - 1]], h, f)[0]); + h_accum += h; + h_list.push(h_accum); + } - for (let i = 0; i < h_list.length; i++) { - result.exact.push(g(h_list[i])); - } + for (let i = 0; i < h_list.length; i++) { + exact.push(g(h_list[i])[0]); + } - result.error = calculate_error(result.estimate, result.exact); + error = calculate_error(estimate, exact); - write_table(result.estimate, result.exact, result.error, h_list, name, 10, error_list); + write_table(estimate, exact, error, h_list, name, 10, error_list); + if (adaptive) { 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); - + } else { + error.forEach(e => { + expect(e).toBeLessThan(1); + }); } +} +describe("One-dimensional integration", () => { + // exact solution + function g(x: number): Array<number> { + return [2 * math.exp(-0.5 * x)]; + }; -// ------------------------------------------------------------------------------------------------ + // differential equation + function f(x: Array<number>): Array<number> { + return [-0.5 * x[0]]; + } -describe("Euler", () => { + describe("Euler", () => { - it("Fixed step size", () => { - fixed_step_size(Integrator.Euler, 0.1, f, "Euler-Fixed"); - }); + it("Fixed step size", () => { + eval_1d_integrator(Integrator.Euler, 0.1, f, g, "Euler-Fixed"); + }); - it("Adaptive step size", () => { - adaptive_step_size(Integrator.Euler, 0.1, f, "Euler-Adaptive"); - }); -}) + it("Adaptive step size", () => { + eval_1d_integrator(Integrator.Euler, 0.1, f, g, "Euler-Adaptive", true); + }); + }) -describe("Midpoint", () => { - it("Fixed step size", () => { - fixed_step_size(Integrator.Midpoint, 0.1, f, "Midpoint-Fixed"); - }); + describe("Midpoint", () => { + it("Fixed step size", () => { + eval_1d_integrator(Integrator.Midpoint, 0.1, f, g, "Midpoint-Fixed"); + }); - it("Adaptive step size", () => { - adaptive_step_size(Integrator.Midpoint, 0.1, f, "Midpoint-Adaptive"); + it("Adaptive step size", () => { + eval_1d_integrator(Integrator.Midpoint, 0.1, f, g, "Midpoint-Adaptive", true); + }); }); + + describe("Runge-Kutta", () => { + it("Fixed step size", () => { + eval_1d_integrator(Integrator.Runge_Kutta, 0.1, f, g, "Runge-Kutta-Fixed"); + }); + + it("Adaptive step size", () => { + eval_1d_integrator(Integrator.Runge_Kutta, 0.1, f, g, "Runge-Kutta-Adaptive", true); + }); + }) }); -describe("Runge-Kutta", () => { - it("Fixed step size", () => { - fixed_step_size(Integrator.Runge_Kutta, 0.1, f, "Runge-Kutta-Fixed"); +// ------------------------------------------------------------------------------------------------ + +function eval_2d_integrator(integrator: (x: Array<number>, h: number, f: (x: Array<number>) => Array<number>) => Array<number>, +h: number, f: (x: Array<number>) => Array<number>, g: (x: number) => Array<number>, name: string, adaptive: boolean = false): void { + + let h_list: Array<number> = []; + let error_list: Array<number> = []; + + let error: Array<Array<number>> = []; + let exact: Array<Array<number>> = []; + let estimate: Array<Array<number>> = []; + + h_list.push(0); + let h_accum = h; + estimate.push(g(h_list[0])); + while (h_accum < x_max) { + if (adaptive) { + h = Integrator.Adaptive_step_size(integrator, estimate[estimate.length - 1], h, f, error_max, (adaptive) ? error_list : undefined); + } + estimate.push(integrator(estimate[estimate.length - 1], h, f)); + h_accum += h; + h_list.push(h_accum); + } + + for (let i = 0; i < h_list.length; i++) { + exact.push(g(h_list[i])); + } + + error = calculate_error_2d(estimate, exact); + + write_table_2d(estimate, exact, error, h_list, name, 5, error_list); + + +} + +describe("Mulit-dimensional integration", () => { + + const a = [0, 9.81]; + const p0 = [0, 0]; + const v0 = [1, 1]; + + // exact solution + function g(t: number): Array<number> { + return [a[0] / 2 * t * t + v0[0] * t + p0[0], a[1] / 2 * t * t + v0[1] * t + p0[1], a[0] * t + v0[0], a[1] * t + v0[1]]; + } + + // differential equation + function f(x: Array<number>): Array<number> { + return [x[2], x[3], a[0], a[1]]; + } + + describe("Euler", () => { + it("Fixed step size", () => { + eval_2d_integrator(Integrator.Euler, 0.1, f, g, "2D-Euler-Fixed"); + }); + + it("Adaptive step size", () => { + eval_2d_integrator(Integrator.Euler, 0.1, f, g, "2D-Euler-Adaptive", true); + }); + }) + + describe("Midpoint", () => { + it("Fixed step size", () => { + eval_2d_integrator(Integrator.Midpoint, 0.1, f, g, "2D-Midpoint-Fixed"); + }); + + it("Adaptive step size", () => { + eval_2d_integrator(Integrator.Midpoint, 0.1, f, g, "2D-Midpoint-Adaptive", true); + }); }); - it("Adaptive step size", () => { - adaptive_step_size(Integrator.Runge_Kutta, 0.1, f, "Runge-Kutta-Adaptive"); + describe("Runge-Kutta", () => { + it("Fixed step size", () => { + eval_2d_integrator(Integrator.Runge_Kutta, 0.1, f, g, "2D-Runge-Kutta-Fixed"); + }); + + it("Adaptive step size", () => { + eval_2d_integrator(Integrator.Runge_Kutta, 0.1, f, g, "2D-Runge-Kutta-Adaptive", true); + }); }); -}) \ No newline at end of file +});