JS 中的龙格库塔问题
Runge Kutta problems in JS
我正在尝试对 Javascript 中 spring 的质量进行 Runge-Kutta 实现,并使用 D3 将其可视化。目的是将其与正向欧拉进行比较并评论差异。我的 FE 工作正常,绘图也很好,但 Runge-Kutta 朝负方向发射并且从不环绕。
这是一个 plunkr 的 vis 和代码,但我也会添加 JS(仅用于 ODE 求解器)。
// *** Functions for ODE Solvers *** //
function FEx (x, v, h)
{
return x + h*v;
}
function FEv (x, v, h)
{
var k = 1; var m = 0.5; var g = 0;
return v + h*( (-k/m)*x + g );
}
function RKx (x, v, h)
{
var k1 = FEx(x, v, h);
var k2 = FEx(x+h/2*k1, v+h/2, h);
var k3 = FEx(x+h/2*k2, v+h/2, h);
var k4 = FEx(x+h*k3, v+h, h);
return x + h/6*(k1 + 2*k2 + 2*k3 + k4);
}
function RKy (x, v, h)
{
var k1 = FEv(x, v, h);
var k2 = FEv(x+h/2, v+h/2*k1, h);
var k3 = FEv(x+h/2, v+h/2*k2, h);
var k4 = FEv(x+h, v+h*k3, h);
return v + h/6*(k1 + 2*k2 + 2*k3 + k4);
}
// FORWARD EULER
function forewardEuler (x, v, h, n)
{
// Initialize an array to hold the values
// JS doesn't really support multi-dimensional arrays
// so this is a "jagged" nested array
var values = new Array(n);
for(i = 0; i < values.length; i++)
values[i] = new Array(2);
// Initial conditions
values[0] = [x, v];
for (i = 1; i < n; ++i)
{
values[i][0] = FEx(values[i-1][0], values[i-1][1], h);
values[i][1] = FEv(values[i-1][0], values[i-1][1], h);
}
return values;
}
// 4TH ORDER RUNGE-KUTTA
function RK4 (x, v, h, n)
{
// Initialize an array to hold the values
var values = new Array(n);
for(i = 0; i < values.length; i++)
values[i] = new Array(2);
// Initial conditions
values[0] = [x, v];
for (i = 1; i < n; ++i)
{
values[i][0] = RKx(values[i-1][0], values[i-1][1], h);
values[i][1] = RKy(values[i-1][0], values[i-1][1], h);
}
return values;
}
// *** Setting up the data *** //
var rkValues = RK4(1, 0, 0.1, 100);
var feValues = forewardEuler(1, 0, 0.1, 100);
这有一些非常基本的概念问题。对于耦合系统,您必须同时评估所有操作。也就是说,在 y'(t)=f(y(t))
中,函数 y(t)
是矢量值,f
具有矢量作为输入和矢量作为输出。欧拉法则可以概括为
k = f(y[i]);
y[i+1] = y[i] + h*k;
允许灵活地评估 f
的组件。 RK4 遵循类似的方案,斜率 k0,...,k3
是函数 f
在不同修改点的所有值。
欧拉步当然不是 RK4 步的一部分,也不应与 ODE 的系统函数混淆。
所以你应该在
的方向上使用一些东西
function odefuncX(x,v) {return v;}
function odefuncV(x,v) {
var k = 1; var m = 0.5; var g = 0;
return (-k/m)*x + g;
}
function EulerStep(x,v,h) {
var kx = odefuncX(x,v);
var kv = odefuncV(x,v);
return [ x+h*kx, v+h*kv ];
}
function RK4Step(x,v,h) {
var kx0 = odefuncX(x,v);
var kv0 = odefuncV(x,v);
var kx1 = odefuncX(x+0.5*h*kx0,v+0.5*h*kv0);
var kv1 = odefuncV(x+0.5*h*kx0,v+0.5*h*kv0);
var kx2 = odefuncX(x+0.5*h*kx1,v+0.5*h*kv1);
var kv2 = odefuncV(x+0.5*h*kx1,v+0.5*h*kv1);
var kx3 = odefuncX(x+ h*kx2,v+ h*kv2);
var kv3 = odefuncV(x+ h*kx2,v+ h*kv2);
return [ x+h/6*(kx0+2*(kx1+kx2)+kx3),
v+h/6*(kv0+2*(kv1+kv2)+kv3) ];
}
// 4TH ORDER RUNGE-KUTTA
function RK4 (x, v, h, n) {
// Initialize an array to hold the values
var values = new Array(n);
// Initial conditions
values[0] = [x, v];
for (i = 1; i < n; ++i) {
values[i] = RK4Step(values[i-1][0], values[i-1][1], h);
}
return values;
}
我正在尝试对 Javascript 中 spring 的质量进行 Runge-Kutta 实现,并使用 D3 将其可视化。目的是将其与正向欧拉进行比较并评论差异。我的 FE 工作正常,绘图也很好,但 Runge-Kutta 朝负方向发射并且从不环绕。
这是一个 plunkr 的 vis 和代码,但我也会添加 JS(仅用于 ODE 求解器)。
// *** Functions for ODE Solvers *** //
function FEx (x, v, h)
{
return x + h*v;
}
function FEv (x, v, h)
{
var k = 1; var m = 0.5; var g = 0;
return v + h*( (-k/m)*x + g );
}
function RKx (x, v, h)
{
var k1 = FEx(x, v, h);
var k2 = FEx(x+h/2*k1, v+h/2, h);
var k3 = FEx(x+h/2*k2, v+h/2, h);
var k4 = FEx(x+h*k3, v+h, h);
return x + h/6*(k1 + 2*k2 + 2*k3 + k4);
}
function RKy (x, v, h)
{
var k1 = FEv(x, v, h);
var k2 = FEv(x+h/2, v+h/2*k1, h);
var k3 = FEv(x+h/2, v+h/2*k2, h);
var k4 = FEv(x+h, v+h*k3, h);
return v + h/6*(k1 + 2*k2 + 2*k3 + k4);
}
// FORWARD EULER
function forewardEuler (x, v, h, n)
{
// Initialize an array to hold the values
// JS doesn't really support multi-dimensional arrays
// so this is a "jagged" nested array
var values = new Array(n);
for(i = 0; i < values.length; i++)
values[i] = new Array(2);
// Initial conditions
values[0] = [x, v];
for (i = 1; i < n; ++i)
{
values[i][0] = FEx(values[i-1][0], values[i-1][1], h);
values[i][1] = FEv(values[i-1][0], values[i-1][1], h);
}
return values;
}
// 4TH ORDER RUNGE-KUTTA
function RK4 (x, v, h, n)
{
// Initialize an array to hold the values
var values = new Array(n);
for(i = 0; i < values.length; i++)
values[i] = new Array(2);
// Initial conditions
values[0] = [x, v];
for (i = 1; i < n; ++i)
{
values[i][0] = RKx(values[i-1][0], values[i-1][1], h);
values[i][1] = RKy(values[i-1][0], values[i-1][1], h);
}
return values;
}
// *** Setting up the data *** //
var rkValues = RK4(1, 0, 0.1, 100);
var feValues = forewardEuler(1, 0, 0.1, 100);
这有一些非常基本的概念问题。对于耦合系统,您必须同时评估所有操作。也就是说,在 y'(t)=f(y(t))
中,函数 y(t)
是矢量值,f
具有矢量作为输入和矢量作为输出。欧拉法则可以概括为
k = f(y[i]);
y[i+1] = y[i] + h*k;
允许灵活地评估 f
的组件。 RK4 遵循类似的方案,斜率 k0,...,k3
是函数 f
在不同修改点的所有值。
欧拉步当然不是 RK4 步的一部分,也不应与 ODE 的系统函数混淆。
所以你应该在
的方向上使用一些东西function odefuncX(x,v) {return v;}
function odefuncV(x,v) {
var k = 1; var m = 0.5; var g = 0;
return (-k/m)*x + g;
}
function EulerStep(x,v,h) {
var kx = odefuncX(x,v);
var kv = odefuncV(x,v);
return [ x+h*kx, v+h*kv ];
}
function RK4Step(x,v,h) {
var kx0 = odefuncX(x,v);
var kv0 = odefuncV(x,v);
var kx1 = odefuncX(x+0.5*h*kx0,v+0.5*h*kv0);
var kv1 = odefuncV(x+0.5*h*kx0,v+0.5*h*kv0);
var kx2 = odefuncX(x+0.5*h*kx1,v+0.5*h*kv1);
var kv2 = odefuncV(x+0.5*h*kx1,v+0.5*h*kv1);
var kx3 = odefuncX(x+ h*kx2,v+ h*kv2);
var kv3 = odefuncV(x+ h*kx2,v+ h*kv2);
return [ x+h/6*(kx0+2*(kx1+kx2)+kx3),
v+h/6*(kv0+2*(kv1+kv2)+kv3) ];
}
// 4TH ORDER RUNGE-KUTTA
function RK4 (x, v, h, n) {
// Initialize an array to hold the values
var values = new Array(n);
// Initial conditions
values[0] = [x, v];
for (i = 1; i < n; ++i) {
values[i] = RK4Step(values[i-1][0], values[i-1][1], h);
}
return values;
}