Сравнение методов Рунге-Кутта и Липфрога
Материал из Department of Theoretical and Applied Mechanics
Версия от 18:31, 14 января 2016; 31.134.191.119 (обсуждение)
Виртуальная лаборатория>Сравнение методов Рунге-Кутта и leapfrog
Содержание
Постановка задачи
Дано простейшее уравнение движения грузика на пружине: . Необходимо интегрировать его с помощью двух методов: leapfrog и Рунге-Кутта. Построить фазовую плоскость для каждого из методов и сравнить результаты.
Метод численного интегрирования Рунге-Кутты
Для нашего уравнения алгоритм будет выглядеть следующим образом:
Вычисление нового значения проходит в четыре стадии:
где - шаг интегрирования.
Тогда приближенное значение в последующих точках вычисляется по итерационной формуле:
Метод численного интегрирования leapfrog
Для уравнения второй степени скорость и перемещение находятся следующим образом:
Метод Рунге-Кутта Метод leapfrog
Скачать программу Particles1.rar.
Текст программы на языке JavaScript (разработчик Погодина Валерия):
Файл "Particles1.js"
1 <!DOCTYPE html>
2 <html>
3 <head>
4 <meta charset="UTF-8" />
5 <title> Particle</title>
6 <script src="Particles2.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 }
Выводы
Мы наблюдаем симплектический характер метода Липфрога: сохраняет энергию(слегка измененную). Метод Рунге-Кутты напротив не соханяет энергию системы.
Ссылки
- leapfrog
- Рунге-Кутты
- Курсовые работы по ВМДС: 2015-2016
- Введение в механику дискретных сред
- Виртуальная лаборатория
