Сравнение методов Рунге-Кутта и Липфрога — различия между версиями

Материал из Department of Theoretical and Applied Mechanics
Перейти к: навигация, поиск
Строка 6: Строка 6:
  
 
== '''Метод численного интегрирования Рунге-Кутты''' ==
 
== '''Метод численного интегрирования Рунге-Кутты''' ==
−
Рассмотри задачу Коши для системы обыкновенных дифференциальных уравнений первого порядка.
+
Рассмотри задачу Коши для системы обыкновенных дифференциальных уравнений первого порядка.<br />
−
<math>\dot{y} = f(x, y), y(x_0)=y_0</math>
+
<math>\dot{y} = f(x, y), y(x_0)=y_0</math><br />
−
Тогда приближенное значение в последующих точках вычисляется по итерационной формуле:
+
Тогда приближенное значение в последующих точках вычисляется по итерационной формуле:<br />
−
<math>y_n+1 = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)</math>
+
<math>y_n+1 = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)</math><br />
−
=='''Основные уравнения'''==
+
Вычисление нового зачения проходит в четыре стадии:<br />
−
<div style="font-size: 16px;">1)''Расчет статических координат:''</div>
+
<math>k_1 = f(x_n, y_n)</math><br />
−
 
+
<math>k_2 = f(x_n+\frac{h}{2}, y_n+\frac{h}{2}k_1)</math><br />
−
<math>{y}_{i}=i{l}_{0}+\sum^{i}_{k=1} {\Delta l_k}</math>
+
<math>k_3 = f(x_n+\frac{h}{2}, y_n+\frac{h}{2}k_2)</math><br />
−
<br/>
+
<math>k_4 = f(x_n+h, y_n+hk_3)</math><br />
−
<br/>
+
где <math>h</math> - величина сетки по <math>x</math>.<br />  
−
<math>\Delta{l_k} = \frac{(n-k+1)mg}{C} </math>
+
=='''Метод численного интегрирования Липфрога'''==
−
<br/>
+
Для уравнения второй степени скорость и перемещение находятся следующим образом:<br />
−
<br/>
+
<math>\dot{x(t+dt)} = \dot{x(t)}+f(t)dt</math><br />
−
где  <math>{n}</math> — количество грузиков ; <math>{i}</math> — текущий грузик ; <math>{l_0}</math> — длина не растянутой пружины; <math>\Delta{l_k}</math> — статическое удлинение пружины.
+
<math>x(t+dt) = x(t)+\dot{x(t+dt)}dt</math><br />
−
 
+
                            '''Метод Рунге-Кутта'''                                                            '''Метод Липфрога'''
−
<div style="font-size: 16px;">2) ''Уравнение движения грузиков при отсутствии верхней пружины:''</div>
+
{{#widget:Iframe |url=http://tm.spbstu.ru/htmlets/Pogodina/Integrate/particles1.html |width=1200|height=800 |border=0 }}
−
<math>\ddot y_i = - \frac{C}{m} ({y}_{i+1} - 2{y}_{i}+{y}_{i-1}) + g </math> 
 
−
<br/>
 
−
<br/>
 
−
при <math>i = n </math> :
 
−
<br/>
 
−
<math> \ddot y_n = - \frac{C}{m}(y_n - {y}_{n-1} - l_0) + g</math>
 
−
<br/>
 
−
<br/>'''Программа (лучше всего смотреть ее в Mozilla Firefox)'''
 
−
<br/>
 
−
{{#widget:Iframe |url=http://tm.spbstu.ru/htmlets/Alexandrov/Spring_illusion.html |width=800|height=800 |border=0 }}
 
 
<br />
 
<br />
 
<br />
 
<br />
−
Скачать программу [[Медиа:Spring_illusion.rar|Spring_illusion.rar]].
+
Скачать программу [[Медиа:Particles1.rar|Particles1.rar]].
 
<br />
 
<br />
 
<br />
 
<br />
 
<div class="mw-collapsible mw-collapsed" style="width:100%" >
 
<div class="mw-collapsible mw-collapsed" style="width:100%" >
−
'''Текст программы на языке JavaScript (разработчик [[Александров Сергей]]):''' <div class="mw-collapsible-content">
+
'''Текст программы на языке JavaScript (разработчик [[Погодина Валерия]]):''' <div class="mw-collapsible-content">
−
Файл '''"Spring_illusin.js"'''
+
Файл '''"Particles1.js"'''
 
<syntaxhighlight lang="javascript" line start="1" enclose="div">
 
<syntaxhighlight lang="javascript" line start="1" enclose="div">
−
window.addEventListener("load", Main_Spring, true);
+
<!DOCTYPE html>
−
function Main_Spring() {
+
<html>
 +
<head>
 +
    <meta charset="UTF-8" />
 +
    <title> Particle</title>
 +
    <script src="Particles1.js"></script>
 +
<!--    <script src="jquery.min.js"></script>
 +
    <script src="jquery.flot.js"></script> -->
 +
</head>
 +
<body>
 +
<br>
 +
<table><tr valign="top">
 +
<td>
 +
    <canvas id="canvasGraph" width="600" height="600" style="border:1px solid #000000;"></canvas>
 +
<canvas id="canvasGraph1" width="600" height="600" style="border:1px solid #000000;"></canvas>
 +
</td>
 +
<!-- <td>
 +
<input type="button" id="drop_spring" value="Drop Spring"/><br>
 +
<br/>
 +
<input type="button" id="new_static" value="New Static"/><br>
 +
<br>
 +
    <input type="range" id="slider_m" min="0.1" max="1.5" step=0.1 style="width: 150px;" />
 +
    Масса грузиков = <input type="number" id="number_m" min="0.1" max="1.5" step=0.1 style="width: 50px;" /><br>
 +
 +
<input type="range" id="slider_spf" min="1" max="15" step=1 style="width: 150px;" />
 +
    Скорость протекания процесса = <input type="number" id="number_spf" min="1" max="15" step=1 style="width: 50px;" /><br> -->
 +
 +
<input type="range" id="slider_n" min="1" max="100" step=1 style="width: 150px;" />
 +
    Количество частиц = <input type="number" id="number_n" min="1" max="100" step=1 style="width: 50px;" /><br>
 +
<br>
 +
<input type="range" id="slider_dt" min="0.0001" max="0.01" step=0.0001 style="width: 150px;" />
 +
    dt = <input type="number" id="number_dt" min="0.0001" max="0.01" step=0.0001 style="width: 50px;" /><br>
 +
<br>
 +
<script
 +
type="text/javascript"> app = new MainParticle(document.getElementById('canvasGraph'),document.getElementById('canvasGraph1'));
 +
</script>
 +
<!-- График перемещений последнего грузика:
 +
    <div id="vGraph1" style="width:400px; height:200px; clear:both;"></div> -->
 +
  
−
// *** Некие исходные данные ***
+
</tr></table>
 +
</body>
 +
</html>
 +
function MainParticle(canvas_gr, canvas_gr1) {
 +
    // Предварительные установки
 +
var context_gr  = canvas_gr.getContext("2d");  // на context происходит рисование
 +
var context_gr1  = canvas_gr1.getContext("2d");  // на context происходит рисование
 +
  //  Задание констант
 +
    const Pi = 3.1415926;                  // число "пи"
 +
    const m0 = 1;                          // масштаб массы
 +
    const T0 = 1;                          // масштаб времени (период колебаний исходной системы)
 +
    const a0 = 1;                          // масштаб расстояния (диаметр шара)
  
−
var canvas = spring_canvas;
+
    const g0 = a0 / T0 / T0;               // масштаб ускорения (ускорение, при котором за T0 будет пройдено расстояние a0)
−
canvas.onselectstart = function () {return false;}; // запрет выделения canvas
+
    const k0 = 2 * Pi / T0;                 // масштаб частоты
−
var ctx = canvas.getContext("2d"); // на ctx происходит рисование
+
    const C0 = m0 * k0 * k0;               // масштаб жесткости
−
var w = canvas.width; // ширина окна в расчетных координатах
+
    const B0 = 2 * m0 * k0;                 // масштаб вязкости
−
var h = canvas.height; // высота окна в расчетных координатах
 
  
−
var Pi = 3.1415926; // число "пи"
+
    // *** Задание физических параметров ***
−
var g = 9.8; // гравитационная постоянная
 
  
−
var T0 = 0.01; // масштаб времени (период колебаний исходной системы)
+
    const Ny = 5;                          // число шаров, помещающихся по вертикали в окно (задает размер шара относительно размера окна)
 +
    const m = 1 * m0;                      // масса
 +
    const Cwall = 10 * C0;                  // жесткость стен
 +
    const  B = 0.01 * B0;                  // вязкость среды
 +
    const  B1 = 0.03 * B0;                  // вязкость на стенках
 +
    var   mg = 0.25 * m * g0;              // сила тяжести
 +
    const  r = 0.5 * a0;                    // радиус частицы в расчетных координатах
 +
var stiff = 1 * C0;                     // "жесткость" пружинки
  
−
var m0 = 0.15; // масштаб массы маятника
+
var vx0 = 1 * a0 / T0;
−
var C0 = 1; // масштаб жесткости пружины
 
  
−
var count = true; // проводить ли расчет системы
+
−
var v = 0; // скорость тела
+
−
var t = 0;
+
 +
 +
    // *** Задание вычислительных параметров ***
  
−
// параметры полученные из размеров холста
+
    const fps = 50;                         // frames per second - число кадров в секунду (качеcтво отображения)
−
var rw = canvas.width / 4;
+
    const spf = 100;                       // steps per frame  - число шагов интегрирования между кадрами (скорость расчета)
−
var rh = canvas.height / 100;
+
    var dt  = 0.2 * T0 / fps;  
−
 
+
var n = 100; // шаг интегрирования
−
// параметры грузиков
 
−
var x0 = 0;
 
−
var y0 = 0;
 
−
var vy0 = 0;
 
−
var rad0 = 7;
 
−
var sTime = 0;
 
−
// параметры пружины
 
−
var Lp = 30; //длина пружины
 
−
 
 
−
 
 
−
// *** Задание вычислительных параметров ***
 
−
 
 
−
var fps = 60; // frames per second - число кадров в секунду (качечтво отображения)
 
−
var spf = 5; // steps per frame  - число шагов интегрирования между кадрами
 
−
var dt = 50 * T0 / fps; // шаг интегрирования (качество расчета)
 
−
var steps = 0; // количество шагов интегрирования
 
−
 
 
−
 
 
−
// *** Задание физических параметров ***
 
−
var m = 1 * m0; // масса грузика
 
−
var C = 1 * C0; // жесткость пружины
 
−
var L0 = 0; // изначальное удлинение пружины
 
−
var n = 5; // количество тел
 
−
var circle; // круглые тела
 
−
 
 
−
// *** Установка слайдеров для переменных величин ***
 
−
slider_m.value = (m / m0).toFixed(1);
 
−
number_m.value = (m / m0).toFixed(1);
 
−
slider_spf.value = (spf).toFixed(1);
 
−
number_spf.value = (spf).toFixed(1);
 
 
slider_n.value = parseInt(n);
 
slider_n.value = parseInt(n);
 
number_n.value = parseInt(n);
 
number_n.value = parseInt(n);
−
 
+
slider_dt.value = dt;
−
function setM(new_m) {
+
number_dt.value = dt;
−
m = new_m * m0;
+
 +
        function setN(new_n) {
 +
                n = new_n;
 +
for (var i = 0; i< n; i++) {
 +
for (var j = 0; j < n; j++) {
 +
c = [];
 +
d = [];
 +
c.x = w / 2 + w / 7.1 / n * i;
 +
c.vx = vx0 / 1.7 / n * j;
 +
d.x = 0.28 / n * i;
 +
d.vx = vx0 / 1.7 / n * j;
 +
b[i][j] = c;
 +
a[i][j] = d;
 +
}
 
}
 
}
−
 
+
context_gr.clearRect(0, 0, w * scale, h * scale);
−
function setSpf(new_spf) {
+
context_gr1.clearRect(0, 0, w * scale, h * scale);
−
spf = new_spf;
+
        }
 +
 +
function setdt(new_dt) {
 +
dt = new_dt;
 +
for (var i = 0; i< n; i++) {
 +
for (var j = 0; j < n; j++) {
 +
c = [];
 +
d = [];
 +
c.x = w / 2 + w / 7.1 / n * i;
 +
c.vx = vx0 / 1.7 / n * j;
 +
d.x = 0.28 / n * i;
 +
d.vx = vx0 / 1.7 / n * j;
 +
b[i][j] = c;
 +
a[i][j] = d;
 +
}
 
}
 
}
−
function setN(new_n) {
+
context_gr.clearRect(0, 0, w * scale, h * scale);
−
n = new_n;
+
context_gr1.clearRect(0, 0, w * scale, h * scale);
−
StaticStart();
+
        }
−
}
+
−
 
 
−
slider_m.oninput = function () {
 
−
number_m.value = slider_m.value;
 
−
setM(slider_m.value);
 
−
};
 
−
number_m.oninput = function () {
 
−
slider_m.value = number_m.value;
 
−
setM(number_m.value);
 
−
};
 
−
 
 
−
 
 
−
slider_spf.oninput = function () {
 
−
number_spf.value = slider_spf.value;
 
−
setSpf(slider_spf.value);
 
−
};
 
−
number_spf.oninput = function () {
 
−
slider_spf.value = number_spf.value;
 
−
setSpf(number_spf.value);
 
−
};
 
−
 
 
 
slider_n.oninput = function () {
 
slider_n.oninput = function () {
 
number_n.value = slider_n.value;
 
number_n.value = slider_n.value;
 
setN(slider_n.value);
 
setN(slider_n.value);
−
};
+
        };
 +
 
number_n.oninput = function () {
 
number_n.oninput = function () {
 
slider_n.value = number_n.value;
 
slider_n.value = number_n.value;
 
setN(number_n.value);
 
setN(number_n.value);
−
};
+
        };
  
−
// *** Условие падения тела ***
+
slider_dt.oninput = function () {
−
var spring_dropped = false;
+
number_dt.value = slider_dt.value;
−
drop_spring.onclick = function () {
+
setdt(slider_dt.value);
−
spring_dropped = true;
+
        };
−
};
+
 +
number_dt.oninput = function () {
 +
slider_dt.value = number_dt.value;
 +
setdt(number_dt.value);
 +
        };
  
−
// *** Пересчет статических координат
+
−
new_static.onclick = function () {
+
−
StaticStart()
+
// Задание констант для рисования
−
};
+
const scale = canvas_gr.height / Ny / a0;  // масштабный коэффициент для перехода от расчетных к экранным координатам
 +
 +
var w = canvas_gr.width / scale;          // ширина окна в расчетных координатах
 +
    var h = canvas_gr.height / scale;          // высота окна в расчетных координатах
 +
  
−
// *** Создание массива элементов ***
+
    // -------------------------------        Выполнение программы              ------------------------------------------
−
function StaticStart() {
+
// Добавление шара
−
spring_dropped = false;
+
var b = [];
−
circle = [];
+
var a = [];
−
for (var i = 0; i < n; i++) {
+
for (var z = 0; z < n; z++) {
−
var circ = {};
+
b[z] = [];
−
circ.x = x0;
+
a[z] = [];
−
circ.y = y0;
+
}
−
circ.vy = vy0;
+
for (var i = 0; i< n; i++) {
−
circ.rad = rad0;
+
for (var j = 0; j < n; j++) {
−
circ.L = L0;
+
c = [];
−
 
+
d = [];
−
var rgb = HSVtoRGB(i / n * 2, 1, 1);
+
c.x = w / 2 + w / 7.1 / n * i;
−
 
+
c.vx = vx0 / 1.7 / n * j;
−
circ.fill = "rgba(" + rgb.r + ", " + rgb.g + ", " + rgb.b + ", 1)";
+
d.x = 0.28 / n * i;
−
circle[i] = circ;
+
d.vx = vx0 / 1.7 / n * j;
−
}
+
b[i][j] = c;
−
for (var i = 0; i < n; i++) {
+
a[i][j] = d;
−
if (i == 0) {
 
−
circle[i].x = w / 2;
 
−
circle[i].y0 = Lp + (n - i) * m * g / C;
 
−
circle[i].y = circle[i].y0;
 
−
 
 
−
} else {
 
−
circle[i].x = w / 2;
 
−
circle[i].y0 = circle[i - 1].y0 + Lp + (n - i) * m * g / C;
 
−
circle[i].y = circle[i].y0;
 
−
console.log(circle[i].y0);
 
−
 
 
−
}
 
 
}
 
}
 
}
 
}
 
 
−
// *** функция разукрашивания объектов
+
// Основной цикл программы
−
function HSVtoRGB(h, s, v) {             
+
setInterval(control, 1500 / fps);  // функция control вызывается с периодом, определяемым вторым параметром
−
    var r, g, b, i, f, p, q, t;
+
−
    i = Math.floor(h * 6);
+
// ---------------------------------------------------------------------------------------------------------------------
−
    f = h * 6 - i;
+
// ---------------------------------          Определение всех функций              -----------------------------------
−
    p = v * (1 - s);
+
// ---------------------------------------------------------------------------------------------------------------------
−
    q = v * (1 - f * s);
+
−
    t = v * (1 - (1 - f) * s);
+
// основная функция, вызываемая в программе
−
    switch (i % 6) {
+
function control()
−
        case 0: r = v, g = t, b = p; break;
+
{
−
        case 1: r = q, g = v, b = p; break;
+
         physics(); // делаем spf шагов интегрирование
−
        case 2: r = p, g = v, b = t; break;
+
draw_gr(); // рисуем график
−
         case 3: r = p, g = q, b = v; break;
 
−
        case 4: r = t, g = p, b = v; break;
 
−
        case 5: r = v, g = p, b = q; break;
 
 
     }
 
     }
−
    return {
 
−
        r: Math.floor(r * 255),
 
−
        g: Math.floor(g * 255),
 
−
        b: Math.floor(b * 255)
 
−
    };
 
−
}
 
  
−
    // график
+
−
     var vGraph1 = new TM_graph(                 // определить график
+
    // Функция, делающая spf шагов интегрирования
−
        "#vGraph1",                              // на html-элементе #vGraph
+
     function physics() {                    // то, что происходит каждый шаг времен
−
        250,                                    // сколько шагов по оси "x" отображается
+
for (var s = 1; s <= spf; s++) {
−
        0, 100, 5);                           // мин. значение оси Y, макс. значение оси Y, шаг по оси Y
+
for (var k = 0; k < n; k++) {
 +
for (var p = 0; p < n; p++) {
 +
k1v = dt*(-( b[k][p].x - w/2));
 +
k2v = dt*(-( b[k][p].x - w/2)  + k1v/2);
 +
k3v = dt*(-( b[k][p].x - w/2)  + k2v/2);
 +
k4v = dt*(-( b[k][p].x - w/2) + k3v);
 +
 
 +
k1r = dt* b[k][p].vx;
 +
k2r = dt*(b[k][p].vx + k1r/2);
 +
k3r = dt*(b[k][p].vx + k2r/2);
 +
k4r = dt*(b[k][p].vx + k3r);
 +
 
 +
b[k][p].vx += (k1v + 2*k2v + 2*k3v + k4v)/6;
 +
b[k][p].x  += (k1r + 2*k2r + 2*k3r + k4r)/6;
 +
 +
 +
a[k][p].vx = a[k][p].vx - a[k][p].x * dt;         
 +
a[k][p].x = a[k][p].x + a[k][p].vx * dt;
 +
        }
 +
}
 +
}
 +
}
 +
  
−
// *** Функция обеспечивающая "жизнь" пружин ***
+
−
function control() {
+
−
calculate();
 
−
draw();
 
−
requestAnimationFrame(control);
 
−
}
 
−
StaticStart();
 
−
control();
 
  
−
 
+
// Определение функции, рисующей график
−
 
+
context_gr.fillStyle = "#3070d0"; // цвет
−
// *** Функция расчетов координат ***
+
context_gr.strokeStyle = "#ff0000";
−
function calculate() {
+
context_gr1.fillStyle = "#3070d0"; // цвет
−
 
+
context_gr1.strokeStyle = "#ff0000";
−
for (var s = 1; s <= spf; s++) {
+
function draw_gr()  
−
var k = n;
+
{
−
for (var i = 1; i < n; i++) {
+
context_gr.clearRect(0, 0, w * scale, h * scale);
−
if (spring_dropped)
+
context_gr.beginPath();
−
circle[0].L = 0;
+
context_gr.strokeStyle = "#000000";
−
else
+
−
circle[0].L = circle[0].y - Lp;
+
context_gr1.clearRect(0, 0, w * scale, h * scale);
−
circle[i].L = circle[i].y - circle[i - 1].y - Lp;
+
context_gr1.beginPath();
−
}
+
context_gr1.strokeStyle = "#000000";
−
 
+
−
for (var i = 0; i < n; i++) {
+
// ось
−
if (i == n - 1)
+
context_gr.moveTo(w/2*scale, (h)*scale);
−
circle[i].vy += ((-1) * C * (circle[i].L) / m + g) * dt;
+
context_gr.lineTo(w/2*scale, -h*scale);
−
else
+
context_gr1.moveTo(w/2*scale, (h)*scale);
−
circle[i].vy += ((-1) * C * (circle[i].L - circle[i + 1].L) / m + g) * dt;
+
context_gr1.lineTo(w/2*scale, -h*scale);
 +
 +
// ось
 +
context_gr.moveTo(0, (h/2)*scale);
 +
context_gr.lineTo((w)*scale, (h/2)*scale);
 +
context_gr1.moveTo(0, (h/2)*scale);
 +
context_gr1.lineTo((w)*scale, (h/2)*scale);
 +
 +
        context_gr.closePath();
 +
context_gr.stroke();
 +
        context_gr1.closePath();
 +
context_gr1.stroke();
 +
 +
// график
 +
for (var l = 0; l< n; l++) {
 +
for (var m = 0; m < n; m++) {
 +
context_gr.beginPath();
 +
context_gr.arc(b[l][m].x * scale, (-b[l][m].vx + h / 2) * scale, 1, 0, 2 * Math.PI, false);                  //рисуем шар
 +
context_gr.fill();
 +
context_gr.closePath();
 +
 +
context_gr1.beginPath();
 +
context_gr1.arc(a[l][m].x * 300 + 300, (-a[l][m].vx + h / 2) * scale, 1, 0, 2 * Math.PI, false);                  //рисуем шар
 +
context_gr1.fill();
 +
context_gr1.closePath();
 
 
−
+
}
−
circle[i].y += circle[i].vy * dt;
 
−
 
 
 
}
 
}
−
steps++;
 
−
 
−
if (spring_dropped)
 
−
{ sTime ++;
 
−
if (sTime % 50 == 0) vGraph1.graphIter(sTime/50, (circle[n-1].y-circle[n-1].y0)) ;
 
−
};
 
−
}
 
−
}
 
−
 
−
 
−
// *** Функция рисования объектов ***
 
−
function draw() {
 
−
ctx.clearRect(0, 0, w, h);
 
−
 
−
// Рисование закрепления
 
−
ctx.lineWidth = 6;
 
−
ctx.strokeStyle = "#7394cb";
 
−
ctx.beginPath();
 
−
ctx.moveTo(200, 0);
 
−
ctx.lineTo(300, 0);
 
−
ctx.stroke();
 
−
 
−
for (var i = 0; i < n; i++) {
 
−
 
−
// Рисование пружинок
 
−
if (i == 0) {
 
−
if (!spring_dropped)
 
−
draw_spring(0, circle[i].y, w / 2, 10, 40);
 
−
} else
 
−
draw_spring(circle[i - 1].y, circle[i].y, w / 2, 10, 40);
 
−
 
−
 
−
}
 
−
for (var i = 0; i < n; i++) {
 
 
 
−
// Рисование грузиков
+
}
−
ctx.lineWidth = 6;
 
−
ctx.strokeStyle = "#8B008B";
 
−
ctx.beginPath();
 
−
ctx.moveTo(circle[i].x-20, circle[i].y);
 
−
ctx.lineTo(circle[i].x+20, circle[i].y);
 
−
ctx.stroke();
 
−
}
 
−
 
 
−
}
 
−
 
 
−
// *** Функция рисования пружины ***
 
−
function draw_spring(x_start, x_end, y, n, h) {
 
−
ctx.lineWidth = 1;
 
−
 
 
−
var L = x_end - x_start;
 
−
for (var i = 0; i < n; i++) {
 
−
var rgb = HSVtoRGB(i / n * 0.5, 1, 1);
 
−
ctx.strokeStyle = "rgba(" + rgb.r + ", " + rgb.g + ", " + rgb.b + ", 1)";
 
−
 
 
−
var x_st = x_start + L / n * i;
 
−
var x_end = x_start + L / n * (i + 1);
 
−
var l = x_end - x_st;
 
−
ctx.beginPath();
 
−
ctx.bezierCurveTo(y, x_st, y + h, x_st + l / 4, y, x_st + l / 2);
 
−
ctx.bezierCurveTo(y, x_st + l / 2, y - h, x_st + 3 * l / 4, y, x_st + l);
 
−
ctx.stroke();
 
−
}
 
−
}
 
 
}
 
}
 
</syntaxhighlight>
 
</syntaxhighlight>
Строка 317: Строка 302:
 
</div>
 
</div>
  
−
=='''Физическое объяснение'''==
+
=='''Выводы'''==
−
Мы наблюдаем иллюзию "зависания" нижнего грузика потому, что до него должна дойти волна возмущений верхней пружины, прежде чем он придет в движение.
+
Мы наблюдаем симплектический характер метода Липфрога: сохраняет энергию(слегка измененную). Метод Рунге-Кутты напротив не соханяет энергию системы.
  
 
=='''Ссылки'''==
 
=='''Ссылки'''==

Версия 21:03, 12 января 2016

Виртуальная лаборатория>Сравнение методов Рунге-Кутта и Липфрога

Постановка задачи

Дано простейшее уравнение движения грузика на пружине: [math]\ddot{x} = -cx[/math]. Необходимо интегрировать его с помощью двух методов: Липфрога и Рунге-Кутта. Сравнить результаты.

Метод численного интегрирования Рунге-Кутты

Рассмотри задачу Коши для системы обыкновенных дифференциальных уравнений первого порядка.
[math]\dot{y} = f(x, y), y(x_0)=y_0[/math]
Тогда приближенное значение в последующих точках вычисляется по итерационной формуле:
[math]y_n+1 = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)[/math]
Вычисление нового зачения проходит в четыре стадии:
[math]k_1 = f(x_n, y_n)[/math]
[math]k_2 = f(x_n+\frac{h}{2}, y_n+\frac{h}{2}k_1)[/math]
[math]k_3 = f(x_n+\frac{h}{2}, y_n+\frac{h}{2}k_2)[/math]
[math]k_4 = f(x_n+h, y_n+hk_3)[/math]
где [math]h[/math] - величина сетки по [math]x[/math].

Метод численного интегрирования Липфрога

Для уравнения второй степени скорость и перемещение находятся следующим образом:
[math]\dot{x(t+dt)} = \dot{x(t)}+f(t)dt[/math]
[math]x(t+dt) = x(t)+\dot{x(t+dt)}dt[/math]

                           Метод Рунге-Кутта                                                            Метод Липфрога



Скачать программу Particles1.rar.

Текст программы на языке JavaScript (разработчик Погодина Валерия):

Файл "Particles1.js" 

  1 <!DOCTYPE html>
  2 <html>
  3 <head>
  4     <meta charset="UTF-8" />
  5     <title> Particle</title>
  6     <script src="Particles1.js"></script>
  7 <!--     <script src="jquery.min.js"></script>
  8     <script src="jquery.flot.js"></script> -->
  9 </head>
 10 <body>
 11 <br>
 12 <table><tr valign="top">
 13 	<td>
 14     <canvas id="canvasGraph" width="600" height="600" style="border:1px solid #000000;"></canvas>
 15 	<canvas id="canvasGraph1" width="600" height="600" style="border:1px solid #000000;"></canvas>
 16 	</td>
 17 <!-- 	<td>
 18 	<input type="button" id="drop_spring" value="Drop Spring"/><br>
 19 	<br/>
 20 	<input type="button" id="new_static" value="New Static"/><br>
 21 	<br>
 22     <input type="range" id="slider_m" min="0.1" max="1.5" step=0.1 style="width: 150px;" />
 23     Масса грузиков = <input type="number" id="number_m" min="0.1" max="1.5" step=0.1 style="width: 50px;" /><br>
 24 	
 25 		<input type="range" id="slider_spf" min="1" max="15" step=1 style="width: 150px;" />
 26     Скорость протекания процесса = <input type="number" id="number_spf" min="1" max="15" step=1 style="width: 50px;" /><br> -->
 27 	
 28 	<input type="range" id="slider_n" min="1" max="100" step=1 style="width: 150px;" />
 29     Количество частиц = <input type="number" id="number_n" min="1" max="100" step=1 style="width: 50px;" /><br>
 30 	<br>
 31 	<input type="range" id="slider_dt" min="0.0001" max="0.01" step=0.0001 style="width: 150px;" />
 32     dt = <input type="number" id="number_dt" min="0.0001" max="0.01" step=0.0001 style="width: 50px;" /><br>
 33 	<br>
 34 	<script 
 35 	type="text/javascript"> app = new MainParticle(document.getElementById('canvasGraph'),document.getElementById('canvasGraph1'));
 36 	</script>
 37 <!-- 	График перемещений последнего грузика:
 38     <div id="vGraph1" style="width:400px; height:200px; clear:both;"></div> -->
 39 			
 40 
 41 </tr></table>
 42 </body>
 43 </html>
 44 function MainParticle(canvas_gr, canvas_gr1) {
 45     // Предварительные установки
 46 	var context_gr  = canvas_gr.getContext("2d");  // на context происходит рисование
 47 	var context_gr1  = canvas_gr1.getContext("2d");  // на context происходит рисование
 48    //  Задание констант	
 49     const Pi = 3.1415926;                   // число "пи"
 50     const m0 = 1;                           // масштаб массы
 51     const T0 = 1;                           // масштаб времени (период колебаний исходной системы)
 52     const a0 = 1;                           // масштаб расстояния (диаметр шара)
 53 
 54     const g0 = a0 / T0 / T0;                // масштаб ускорения (ускорение, при котором за T0 будет пройдено расстояние a0)
 55     const k0 = 2 * Pi / T0;                 // масштаб частоты
 56     const C0 = m0 * k0 * k0;                // масштаб жесткости
 57     const B0 = 2 * m0 * k0;                 // масштаб вязкости
 58 
 59     // *** Задание физических параметров ***
 60 
 61     const Ny = 5;                           // число шаров, помещающихся по вертикали в окно (задает размер шара относительно размера окна)
 62     const m = 1 * m0;                       // масса
 63     const Cwall = 10 * C0;                  // жесткость стен
 64     const   B = 0.01 * B0;                   // вязкость среды
 65     const  B1 = 0.03 * B0;                   // вязкость на стенках
 66     var    mg = 0.25 * m * g0;               // сила тяжести
 67     const   r = 0.5 * a0;                     // радиус частицы в расчетных координатах
 68 	var stiff = 1 * C0;                     // "жесткость" пружинки
 69 
 70 	var vx0 = 1 * a0 / T0;
 71 
 72 	
 73 	
 74 	
 75 	
 76     // *** Задание вычислительных параметров ***
 77 
 78     const fps = 50;                         // frames per second - число кадров в секунду (качеcтво отображения)
 79     const spf = 100;                        // steps per frame   - число шагов интегрирования между кадрами (скорость расчета)
 80     var dt  = 0.2 * T0 / fps; 
 81 	var n = 100;	// шаг интегрирования 
 82 	slider_n.value = parseInt(n);
 83 	number_n.value = parseInt(n);
 84 	slider_dt.value = dt;
 85 	number_dt.value = dt;
 86 	
 87 	        function setN(new_n) {
 88                 n = new_n;
 89 	for (var i = 0; i< n; i++) {
 90 		for (var j = 0; j < n; j++) {
 91 			c = [];
 92 			d = [];
 93 			c.x = w / 2 + w / 7.1 / n * i;
 94 			c.vx = vx0 / 1.7 / n * j;
 95 			d.x = 0.28 / n * i;
 96 			d.vx = vx0 / 1.7 / n * j;
 97 			b[i][j] = c;
 98 			a[i][j] = d;
 99 		}
100 	}
101 		context_gr.clearRect(0, 0, w * scale, h * scale);
102 		context_gr1.clearRect(0, 0, w * scale, h * scale);
103         }
104 		
105 		function setdt(new_dt) {
106 			dt = new_dt;
107 	for (var i = 0; i< n; i++) {
108 		for (var j = 0; j < n; j++) {
109 			c = [];
110 			d = [];
111 			c.x = w / 2 + w / 7.1 / n * i;
112 			c.vx = vx0 / 1.7 / n * j;
113 			d.x = 0.28 / n * i;
114 			d.vx = vx0 / 1.7 / n * j;
115 			b[i][j] = c;
116 			a[i][j] = d;
117 		}
118 	}
119 		context_gr.clearRect(0, 0, w * scale, h * scale);
120 		context_gr1.clearRect(0, 0, w * scale, h * scale);
121         }
122 		
123 	slider_n.oninput = function () {
124 		number_n.value = slider_n.value;
125 		setN(slider_n.value);
126         };
127 		
128 	number_n.oninput = function () {
129 		slider_n.value = number_n.value;
130 		setN(number_n.value);
131         };
132 
133 	slider_dt.oninput = function () {
134 		number_dt.value = slider_dt.value;
135 		setdt(slider_dt.value);
136         };
137 		
138 	number_dt.oninput = function () {
139 		slider_dt.value = number_dt.value;
140 		setdt(number_dt.value);
141         };
142 
143 	
144 	
145 	// Задание констант для рисования
146 	const scale = canvas_gr.height / Ny / a0;  // масштабный коэффициент для перехода от расчетных к экранным координатам
147 	
148 	var w = canvas_gr.width / scale;           // ширина окна в расчетных координатах
149     var h = canvas_gr.height / scale;          // высота окна в расчетных координатах
150 	
151 
152     // -------------------------------         Выполнение программы              ------------------------------------------
153 	// Добавление шара
154 	var b = [];
155 	var a = [];
156 	for (var z = 0; z < n; z++) {
157 		b[z] = [];
158 		a[z] = [];
159 	}
160 	for (var i = 0; i< n; i++) {
161 		for (var j = 0; j < n; j++) {
162 			c = [];
163 			d = [];
164 			c.x = w / 2 + w / 7.1 / n * i;
165 			c.vx = vx0 / 1.7 / n * j;
166 			d.x = 0.28 / n * i;
167 			d.vx = vx0 / 1.7 / n * j;
168 			b[i][j] = c;
169 			a[i][j] = d;
170 		}
171 	}
172 	
173 	// Основной цикл программы
174 	setInterval(control, 1500 / fps);  // функция control вызывается с периодом, определяемым вторым параметром
175 	
176 // ---------------------------------------------------------------------------------------------------------------------
177 // ---------------------------------           Определение всех функций              -----------------------------------
178 // ---------------------------------------------------------------------------------------------------------------------
179 	
180 	// основная функция, вызываемая в программе
181 	function control() 
182 	{
183         physics(); // делаем spf шагов интегрирование
184 		draw_gr(); // рисуем график
185     }
186 
187 	
188     // Функция, делающая spf шагов интегрирования
189     function physics() {                    // то, что происходит каждый шаг времен	
190 		for (var s = 1; s <= spf; s++) {
191 			for (var k = 0; k < n; k++) {
192 				for (var p = 0; p < n; p++) {	
193 			k1v = dt*(-( b[k][p].x - w/2));
194 			k2v = dt*(-( b[k][p].x - w/2)  + k1v/2);
195 			k3v = dt*(-( b[k][p].x - w/2)  + k2v/2);
196 			k4v = dt*(-( b[k][p].x - w/2) + k3v);
197 			   
198 			k1r = dt* b[k][p].vx;
199 			k2r = dt*(b[k][p].vx + k1r/2);
200 			k3r = dt*(b[k][p].vx + k2r/2);
201 			k4r = dt*(b[k][p].vx + k3r);
202 			   
203 			b[k][p].vx += (k1v + 2*k2v + 2*k3v + k4v)/6;
204 			b[k][p].x  += (k1r + 2*k2r + 2*k3r + k4r)/6;
205 			
206 					
207 			a[k][p].vx = a[k][p].vx - a[k][p].x * dt;          
208 			a[k][p].x = a[k][p].x + a[k][p].vx * dt; 
209         }
210 				}
211 			}			
212 		}
213 		
214 
215 	
216 		
217 
218 	// Определение функции, рисующей график
219 	context_gr.fillStyle = "#3070d0"; // цвет
220 	context_gr.strokeStyle = "#ff0000";
221 	context_gr1.fillStyle = "#3070d0"; // цвет
222 	context_gr1.strokeStyle = "#ff0000";
223 	function draw_gr() 
224 	{
225 		context_gr.clearRect(0, 0, w * scale, h * scale);
226 		context_gr.beginPath();
227 		context_gr.strokeStyle = "#000000";
228 		
229 		context_gr1.clearRect(0, 0, w * scale, h * scale);
230 		context_gr1.beginPath();
231 		context_gr1.strokeStyle = "#000000";
232 		
233 		// ось
234 		context_gr.moveTo(w/2*scale, (h)*scale);
235 		context_gr.lineTo(w/2*scale, -h*scale);
236 		context_gr1.moveTo(w/2*scale, (h)*scale);
237 		context_gr1.lineTo(w/2*scale, -h*scale);
238 		
239 		// ось 
240 		context_gr.moveTo(0, (h/2)*scale);
241 		context_gr.lineTo((w)*scale, (h/2)*scale);
242 		context_gr1.moveTo(0, (h/2)*scale);
243 		context_gr1.lineTo((w)*scale, (h/2)*scale);
244 		
245         context_gr.closePath();
246 		context_gr.stroke();
247         context_gr1.closePath();
248 		context_gr1.stroke();
249 		
250 		// график 
251 		for (var l = 0; l< n; l++) {
252 			for (var m = 0; m < n; m++) {
253 				context_gr.beginPath();
254 				context_gr.arc(b[l][m].x * scale, (-b[l][m].vx + h / 2) * scale, 1, 0, 2 * Math.PI, false);                  //рисуем шар
255 				context_gr.fill();
256 				context_gr.closePath();
257 				
258 				context_gr1.beginPath();
259 				context_gr1.arc(a[l][m].x * 300 + 300, (-a[l][m].vx + h / 2) * scale, 1, 0, 2 * Math.PI, false);                  //рисуем шар
260 				context_gr1.fill();
261 				context_gr1.closePath();
262 				
263 			}
264 		}
265 	
266 	}	
267 }

Выводы

Мы наблюдаем симплектический характер метода Липфрога: сохраняет энергию(слегка измененную). Метод Рунге-Кутты напротив не соханяет энергию системы.

Ссылки


Источник — «http://tm.spbstu.ru/?title=Сравнение_методов_Рунге-Кутта_и_Липфрога&oldid=38742»