Сравнение методов численного интегрирования ОДУ (Рунге-Кутта и leapfrog)

Материал из Department of Theoretical and Applied Mechanics
Версия от 12:59, 24 августа 2016; Wikiadmin2 (обсуждение | вклад) (Wikiadmin2 переименовал страницу Сравнение методов численного интегрирования ОДУ(Рунге-Кутта и leapfrog) в [[Сравнение методов численного инте…)

(разн.) ← Предыдущая | Текущая версия (разн.) | Следующая → (разн.)
Перейти к: навигация, поиск
Виртуальная лаборатория>Сравнение методов численного интегрирования ОДУ(Рунге-Кутта и leapfrog)

Постановка задачи [править]

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

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

Для нашего уравнения алгоритм будет выглядеть следующим образом:
Вычисление нового значения проходит в четыре стадии:
[math]kv_1 = dt(-x(t))[/math]
[math]kv_2 = dt(-x(t)+\frac{kv_1}{2})[/math]
[math]kv_3 = dt(-x(t)+\frac{kv_2}{2})[/math]
[math]kv_4 = dt(-x(t)+kv_3)[/math]

[math]kr_1 = dt\dot{x}(t)[/math]
[math]kr_2 = dt(\dot{x}(t)+\frac{kr_1}{2})[/math]
[math]kr_3 = dt(\dot{x}(t)+\frac{kr_2}{2})[/math]
[math]kr_4 = dt(\dot{x}(t)+kv_3)[/math]
где [math]dt[/math] - шаг интегрирования.
Тогда приближенное значение в последующих точках вычисляется по итерационной формуле:
[math]\dot{x}(t+dt) = \dot{x}(t) + \frac{dt}{6}(kr_1 + 2kr_2 + 2kr_3 + kr_4)[/math]
[math]x(t+dt) = x(t) + \frac{dt}{6}(kv_1 + 2kv_2 + 2kv_3 + kv_4)[/math]

Метод численного интегрирования leapfrog[править]

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

                           Метод Рунге-Кутта                                                            Метод 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: метод почти сохраняет энергию. Метод Рунге-Кутты напротив не сохраняет энергию системы.

Ссылки[править]


Источник — «http://tm.spbstu.ru/?title=Сравнение_методов_численного_интегрирования_ОДУ_(Рунге-Кутта_и_leapfrog)&oldid=43570»