JS 中的龙格库塔问题

问题描述 投票：0回答：1

我正在尝试用 Javascript 对弹簧上的质量进行 Runge-Kutta 实现，并使用 D3 对其进行可视化。目的是将其与前向欧拉进行比较并评论差异。我的有限元效果很好，绘图也很好，但龙格-库塔朝负方向发射，并且永远不会环绕。

这是一个带有 vis 和代码的 plunkr，但我也会添加 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);
javascript numerical-methods ode runge-kutta
1个回答
4
投票

这有一些非常基本的概念问题。对于耦合系统，您必须同时评估所有操作。也就是说，在 

y'(t)=f(y(t))
中，函数 
y(t)
是向量值，
f
将向量作为输入，将向量作为输出。欧拉方法可以总结为

k = f(y[i]);
y[i+1] = y[i] + h*k;

允许灵活评估

f
的组成部分。 RK4 遵循类似的方案，斜率 
k0,...,k3
都是函数 
f
在各个修改点的值。

Euler 步骤当然不是 RK4 步骤的一部分，也不应该与 ODE 的系统函数混淆。

所以你应该按照

的方向使用一些东西 
function EulerStep(x,v,h,odefnX,odefnV) {
    var kx = odefnX(x,v);
    var kv = odefnV(x,v);
    return [ x+h*kx, v+h*kv ];
}

function RK4Step(x,v,h,odefnX,odefnV) {
    var kx0 = odefnX(x,v);
    var kv0 = odefnV(x,v);

    var kx1 = odefnX(x+0.5*h*kx0,v+0.5*h*kv0);
    var kv1 = odefnV(x+0.5*h*kx0,v+0.5*h*kv0);

    var kx2 = odefnX(x+0.5*h*kx1,v+0.5*h*kv1);
    var kv2 = odefnV(x+0.5*h*kx1,v+0.5*h*kv1);

    var kx3 = odefnX(x+    h*kx2,v+    h*kv2);
    var kv3 = odefnV(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, odefnX, odefnV) {
    // 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, odefnX, odefnV); 
    }
    return values;
}



// ----------------------------------
// Example usage

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;
}

RK4(0.0, 1.0 , 0.01, 100, odefuncX, odefuncV);

参见 forked Plunker

// *** Chart variables *** //

var margin = {top: 20, right: 20, bottom: 30, left: 50},
        width = 500 - margin.left - margin.right,
        height = 500 - margin.top - margin.bottom;

var x = d3.scale.linear()
    .range([0, width]);

var y = d3.scale.linear()
    .range([height, 0]);

var xAxis = d3.svg.axis()
    .scale(x)
    .orient("bottom");

var yAxis = d3.svg.axis()
    .scale(y)
    .orient("left");

var line = d3.svg.line()
    .x(function(d) { return x(d.X); })
    .y(function(d) { return y(d.V); });

var odeChart = d3.select("#chart-container").append("svg")
    .attr("width", width + margin.left + margin.right)
    .attr("height", height + margin.top + margin.bottom)
  .append("g")
    .attr("transform", "translate(" + margin.left + "," + margin.top + ")");

// *** Functions for ODE Solvers *** //

function odefunc(y) {
  var x=y[0]; var v=y[1];
    var k = 1; var m = 0.5; var g = 0;
    return [v, (-k/m)*x + g];
}

function EulerStep(y,h,odefn) {
    var k = odefn(y);
    return [ y[0]+h*k[0], y[1]+h*k[1] ];
}

function RK4Step(y,h, odefn) {
    var k0 = odefn(y);
 
    var k1 = odefn([y[0]+0.5*h*k0[0],y[1]+0.5*h*k0[1]]);

    var k2 = odefn([y[0]+0.5*h*k1[0],y[1]+0.5*h*k1[1]]);

    var k3 = odefn([y[0]+    h*k2[0],y[1]+    h*k2[1]]);

    return [ y[0]+h/6*(k0[0]+2*(k1[0]+k2[0])+k3[0]),
             y[1]+h/6*(k0[1]+2*(k1[1]+k2[1])+k3[1]) ];
}

// FORWARD EULER
function forwardEuler (x, v, h, n, odefn)
{
    // 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+1);

    // Initial conditions
    values[0] = [x, v];

    for (i = 1; i <= n; ++i)
    {
        values[i] = EulerStep(values[i-1], h,odefn);
    }

    return values;
}

// 4TH ORDER RUNGE-KUTTA 
function RK4 (x, v, h, n, odefn)
{
    // Initialize an array to hold the values
    var values = new Array(n+1);

    // Initial conditions
    values[0] = [x, v];

    for (i = 1; i <= n; ++i)
    {
        values[i] = RK4Step(values[i-1], h, odefn);
    }

    return values;
}

// *** Setting up the data *** //

var rkValues = RK4(1, 0, 0.1, 100, odefunc);
var feValues = forwardEuler(1, 0, 0.025, 400, odefunc);

var rkData = rkValues.map(function(d) {
  return {
     X: +d[0],
     V: +d[1]
  };});

var feData = feValues.map(function(d) {
  return {
     X: +d[0],
     V: +d[1]
  };});

x.domain(d3.extent(feData, function(d) { return d.X; }));
y.domain(d3.extent(feData, function(d) { return d.V; }));

odeChart.append("g")
  .attr("class", "x axis")
  .attr("transform", "translate(0," + height + ")")
  .call(xAxis);

odeChart.append("g")
  .attr("class", "y axis")
  .call(yAxis)

odeChart.append("path")
  .datum(rkData)
  .attr("class", "rkLine")
  .attr("d", line);

odeChart.append("path")
  .datum(feData)
  .attr("class", "feLine")
  .attr("d", line);

body {
    font: 10px sans-serif;
  }

  #chart-container {
    width: 600px;
    margin: 0 auto;
  }

  .axis path,
  .axis line {
    fill: none;
    stroke: #000;
    shape-rendering: crispEdges;
  }

  .x.axis path {
    display: none;
  }

  .rkLine {
    fill: none;
    stroke: steelblue;
    stroke-width: 1.5px;
  }

  .feLine {
    fill: none;
    stroke: red;
    stroke-width: 1.5px;
  }

<script src="https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.17/d3.min.js"></script>
<div id="chart-container"> <h1>Runge Kutta Order 4</h1 ></div>
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.