-
Notifications
You must be signed in to change notification settings - Fork 15
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #1220 from tervay/armdiffeq
- Loading branch information
Showing
13 changed files
with
584 additions
and
881 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,114 @@ | ||
| export type ODEFunction = ( | ||
| t: number, | ||
| y: number[], | ||
| ) => { | ||
| changeRates: number[]; | ||
| shouldStop: boolean; | ||
| }; | ||
|
|
||
| export default class ODESolver { | ||
| constructor( | ||
| private readonly ode: ODEFunction, | ||
| private readonly y0: number[], | ||
| private readonly t0: number, | ||
| private readonly t1: number, | ||
| ) {} | ||
|
|
||
| euler(resolution: number) { | ||
| const h = (this.t1 - this.t0) / resolution; | ||
| const ts = Array.from(Array(resolution + 1), (_, k) => k * h + this.t0); //time series datapoints | ||
| const ys: number[][] = Array.from(Array(resolution + 1), () => | ||
| Array(this.y0.length).fill(0), | ||
| ); | ||
| ys[0] = this.y0; | ||
|
|
||
| for (let i = 0; i < resolution; i++) { | ||
| console.log(i); | ||
| ys[i + 1] = this.ode(ts[i], ys[i]).changeRates.map( | ||
| (x, j) => ys[i][j] + x * h, | ||
| ); //y_n+1 = y_n + dy/dx *h | ||
| } | ||
|
|
||
| return { | ||
| ts: ts, | ||
| ys: ys, | ||
| }; | ||
| } | ||
|
|
||
| rk4(resolution: number) { | ||
| const h = (this.t1 - this.t0) / resolution; | ||
| const ts = Array.from(Array(resolution + 1), (_, k) => k * h + this.t0); //time series datapoints | ||
| const ys = Array.from(Array(resolution + 1), () => | ||
| Array(this.y0.length).fill(0), | ||
| ); | ||
| ys[0] = this.y0; | ||
|
|
||
| if (this.y0.includes(NaN)) | ||
| console.warn("y0 contains invalid starting value", this.y0); | ||
|
|
||
| let stoppingIndex = 0; | ||
| for (let i = 0; i < resolution; i++) { | ||
| stoppingIndex = i; | ||
| let k = this.ode(ts[i], ys[i]); // f(t, y_n) | ||
| const k1 = k.changeRates; | ||
|
|
||
| if (k.shouldStop) { | ||
| break; | ||
| } | ||
|
|
||
| const s1 = ys[i].map((y, j) => y + (k1[j] * h) / 2); | ||
| k = this.ode(ts[i] + h / 2, s1); // f(t + h/2, y_n + k1*h/2) | ||
| const k2 = k.changeRates; | ||
|
|
||
| if (k.shouldStop) { | ||
| break; | ||
| } | ||
|
|
||
| const s2 = ys[i].map((y, j) => y + (k2[j] * h) / 2); | ||
| k = this.ode(ts[i] + h / 2, s2); // f(t + h/2, y_n + k2*h/2) | ||
| const k3 = k.changeRates; | ||
|
|
||
| if (k.shouldStop) { | ||
| break; | ||
| } | ||
|
|
||
| const s3 = ys[i].map((y, j) => y + k3[j] * h); | ||
| k = this.ode(ts[i] + h, s3); // f(t + h, y_n + k3*h) | ||
| const k4 = k.changeRates; | ||
|
|
||
| if (k.shouldStop) { | ||
| break; | ||
| } | ||
|
|
||
| ys[i + 1] = ys[i].map( | ||
| (x, j) => x + (k1[j] / 6 + k2[j] / 3 + k3[j] / 3 + k4[j] / 6) * h, | ||
| ); //y_n+1 = y_n + (k1 +2*k2 + 2*k3 +k4)/6 *h | ||
| } | ||
|
|
||
| return { | ||
| ts: ts.slice(0, stoppingIndex), | ||
| ys: ys.slice(0, stoppingIndex), | ||
| }; | ||
| } | ||
|
|
||
| midpoint(resolution: number) { | ||
| const h = (this.t1 - this.t0) / resolution; | ||
| const ts = Array.from(Array(resolution + 1), (_, k) => k * h + this.t0); //time series datapoints | ||
| const ys = Array.from(Array(resolution + 1), () => | ||
| Array(this.y0.length).fill(0), | ||
| ); | ||
| ys[0] = this.y0; | ||
|
|
||
| for (let i = 0; i < resolution; i++) { | ||
| const k1 = this.ode(ts[i], ys[i]).changeRates; // f(t, y_n) | ||
|
|
||
| const s1 = ys[i].map((y, j) => y + (k1[j] * h) / 2); // y_n + k1 * h/2 | ||
| const k2 = this.ode(ts[i] + h / 2, s1).changeRates; // f(t + h/2, y_n + k1*h/2) | ||
| ys[i + 1] = ys[i].map((x, j) => x + k2[j] * h); //y_n+1 = y_n + k2 *h | ||
| } | ||
| return { | ||
| ts: ts, | ||
| ys: ys, | ||
| }; | ||
| } | ||
| } |
Oops, something went wrong.
9692a5cThere was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Deploy preview for recalc ready!
✅ Preview
https://recalc-cuiqmv4op-tervay.vercel.app
Built with commit 9692a5c.
This pull request is being automatically deployed with vercel-action