Сравнение методов численного интегрирования ОДУ (Рунге-Кутта и leapfrog)
Материал из Department of Theoretical and Applied Mechanics
Версия от 22:23, 16 января 2016; Leranezabudka (обсуждение | вклад)
Виртуальная лаборатория>Сравнение методов численного интегрирования ОДУ(Рунге-Кутта и 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 T0 = 1; // масштаб времени (период колебаний исходной системы)
51 const a0 = 1; // масштаб расстояния (диаметр шара)
52
53 const k0 = 2 * Pi / T0; // масштаб частоты
54
55 // *** Задание физических параметров ***
56
57 const Ny = 5; // число шаров, помещающихся по вертикали в окно (задает размер шара относительно размера окна)
58
59 var vx0 = 1 * a0 / T0;
60
61
62
63
64
65 // *** Задание вычислительных параметров ***
66
67 const fps = 50; // frames per second - число кадров в секунду (качеcтво отображения)
68 const spf = 100; // steps per frame - число шагов интегрирования между кадрами (скорость расчета)
69 var dt = 0.2 * T0 / fps;
70 var n = 40; // шаг интегрирования
71 slider_n.value = parseInt(n);
72 number_n.value = parseInt(n);
73 slider_dt.value = dt;
74 number_dt.value = dt;
75
76 function setN(new_n) {
77 n = new_n;
78 for (var z = 0; z < n; z++) {
79 b[z] = [];
80 a[z] = [];
81 }
82 for (var i = 0; i< n; i++) {
83 for (var j = 0; j < n; j++) {
84 c = [];
85 d = [];
86 c.x = w / 3 / n * i;
87 c.vx = 2 * vx0 / n * j;
88 d.x = w / 3 / n * i;
89 d.vx = 2 * vx0 / n * j;
90 b[i][j] = c;
91 a[i][j] = d;
92 }
93 }
94 context_gr.clearRect(0, 0, w * scale, h * scale);
95 context_gr1.clearRect(0, 0, w * scale, h * scale);
96 }
97
98 function setdt(new_dt) {
99 dt = new_dt;
100 for (var i = 0; i< n; i++) {
101 for (var j = 0; j < n; j++) {
102 c = [];
103 d = [];
104 c.x = w / 3 / n * i;
105 c.vx = 2 * vx0 / n * j;
106 d.x = w / 3 / n * i;
107 d.vx = 2 * vx0 / n * j;
108 b[i][j] = c;
109 a[i][j] = d;
110 }
111 }
112 context_gr.clearRect(0, 0, w * scale, h * scale);
113 context_gr1.clearRect(0, 0, w * scale, h * scale);
114 }
115
116 slider_n.oninput = function () {
117 number_n.value = slider_n.value;
118 setN(slider_n.value);
119 };
120
121 number_n.oninput = function () {
122 slider_n.value = number_n.value;
123 setN(number_n.value);
124 };
125
126 slider_dt.oninput = function () {
127 number_dt.value = slider_dt.value;
128 setdt(slider_dt.value);
129 };
130
131 number_dt.oninput = function () {
132 slider_dt.value = number_dt.value;
133 setdt(number_dt.value);
134 };
135
136
137
138 // Задание констант для рисования
139 const scale = canvas_gr.height / Ny / a0; // масштабный коэффициент для перехода от расчетных к экранным координатам
140
141 var w = canvas_gr.width / scale; // ширина окна в расчетных координатах
142 var h = canvas_gr.height / scale; // высота окна в расчетных координатах
143
144
145 // ------------------------------- Выполнение программы ------------------------------------------
146 // Добавление шара
147 var b = [];
148 var a = [];
149 for (var z = 0; z < n; z++) {
150 b[z] = [];
151 a[z] = [];
152 }
153 for (var i = 0; i< n; i++) {
154 for (var j = 0; j < n; j++) {
155 c = [];
156 d = [];
157 c.x = w / 3 / n * i;
158 c.vx = 2 * vx0 / n * j;
159 d.x = w / 3 / n * i;
160 d.vx = 2 * vx0 / n * j;
161 b[i][j] = c;
162 a[i][j] = d;
163 }
164 }
165
166 // Основной цикл программы
167 setInterval(control, 1500 / fps); // функция control вызывается с периодом, определяемым вторым параметром
168
169 // ---------------------------------------------------------------------------------------------------------------------
170 // --------------------------------- Определение всех функций -----------------------------------
171 // ---------------------------------------------------------------------------------------------------------------------
172
173 // основная функция, вызываемая в программе
174 function control()
175 {
176 physics(); // делаем spf шагов интегрирование
177 draw_gr(); // рисуем график
178 }
179
180
181 // Функция, делающая spf шагов интегрирования
182 function physics() { // то, что происходит каждый шаг времен
183 for (var s = 1; s <= spf; s++) {
184 for (var k = 0; k < n; k++) {
185 for (var p = 0; p < n; p++) {
186
187 k1v = dt*(-( b[k][p].x));
188 k2v = dt*(-( b[k][p].x) + k1v/2);
189 k3v = dt*(-( b[k][p].x) + k2v/2);
190 k4v = dt*(-( b[k][p].x) + k3v);
191
192 k1r = dt* b[k][p].vx;
193 k2r = dt*(b[k][p].vx + k1r/2);
194 k3r = dt*(b[k][p].vx + k2r/2);
195 k4r = dt*(b[k][p].vx + k3r);
196
197 b[k][p].vx += (k1v + 2*k2v + 2*k3v + k4v)/6;
198 b[k][p].x += (k1r + 2*k2r + 2*k3r + k4r)/6;
199
200
201 a[k][p].vx = a[k][p].vx - a[k][p].x * dt;
202 a[k][p].x = a[k][p].x + a[k][p].vx * dt;
203 }
204 }
205 }
206 }
207
208
209
210
211
212 // Определение функции, рисующей график
213 context_gr.fillStyle = "#3070d0"; // цвет
214 context_gr.strokeStyle = "#ff0000";
215 context_gr1.fillStyle = "#3070d0"; // цвет
216 context_gr1.strokeStyle = "#ff0000";
217 function draw_gr()
218 {
219 context_gr.clearRect(0, 0, w * scale, h * scale);
220 context_gr.beginPath();
221 context_gr.strokeStyle = "#000000";
222
223 context_gr1.clearRect(0, 0, w * scale, h * scale);
224 context_gr1.beginPath();
225 context_gr1.strokeStyle = "#000000";
226
227 // ось
228 context_gr.moveTo(w/2*scale, (h)*scale);
229 context_gr.lineTo(w/2*scale, -h*scale);
230 context_gr1.moveTo(w/2*scale, (h)*scale);
231 context_gr1.lineTo(w/2*scale, -h*scale);
232
233 // ось
234 context_gr.moveTo(0, (h/2)*scale);
235 context_gr.lineTo((w)*scale, (h/2)*scale);
236 context_gr1.moveTo(0, (h/2)*scale);
237 context_gr1.lineTo((w)*scale, (h/2)*scale);
238
239 context_gr.closePath();
240 context_gr.stroke();
241 context_gr1.closePath();
242 context_gr1.stroke();
243
244 // график
245 for (var l = 0; l< n; l++) {
246 for (var m = 0; m < n; m++) {
247 context_gr.beginPath();
248 context_gr.arc((b[l][m].x + w/2)* scale, (-b[l][m].vx + h / 2) * scale, 1, 0, 2 * Math.PI, false); //рисуем шар
249 context_gr.fill();
250 context_gr.closePath();
251
252 context_gr1.beginPath();
253 context_gr1.arc((a[l][m].x + w/2)* scale, (-a[l][m].vx + h / 2) * scale, 1, 0, 2 * Math.PI, false); //рисуем шар
254 context_gr1.fill();
255 context_gr1.closePath();
256
257 }
258 }
259
260 }
261 }
Выводы
Мы наблюдаем симплектический характер метода lepfrog: метод почти сохраняет энергию. Метод Рунге-Кутты напротив не сохраняет энергию системы.
Ссылки
- Описание метода leapfrog
- Описание метода Рунге-Кутты
- Курсовые работы по ВМДС: 2015-2016
- Введение в механику дискретных сред
- Виртуальная лаборатория
