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

[[Виртуальная лаборатория]]>[[Сравнение методов численного интегрирования ОДУ(Рунге-Кутта и leapfrog)]] <HR>

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

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

=='''Метод численного интегрирования leapfrog'''==
Для уравнения второй степени скорость и перемещение находятся следующим образом:<br />
<math>\dot{x}(t+dt) = \dot{x}(t)-x(t)dt</math><br />
<math>x(t+dt) = x(t)+\dot{x}(t+dt)dt</math><br />

                            '''Метод Рунге-Кутта'''                                                            '''Метод leapfrog'''
{{#widget:Iframe |url=http://tm.spbstu.ru/htmlets/Pogodina/Integrate/particles2.html |width=1300|height=800 |border=0 }}
<br />
<br />
Скачать программу [[Медиа:Particles1.rar|Particles1.rar]].
<br />
<br />
<div class="mw-collapsible mw-collapsed" style="width:100%" >
'''Текст программы на языке JavaScript (разработчик [[Погодина Валерия]]):''' <div class="mw-collapsible-content">
Файл '''"Particles1.js"'''
<syntaxhighlight lang="javascript" line start="1" enclose="div">
<!DOCTYPE html>
<html>
<head>
    <meta charset="UTF-8" />
    <title> Particle</title>
    <script src="Particles2.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 T0 = 1;                           // масштаб времени (период колебаний исходной системы)
    const a0 = 1;                           // масштаб расстояния (диаметр шара)

    const k0 = 2 * Pi / T0;                 // масштаб частоты

    // *** Задание физических параметров ***

    const Ny = 5;                           // число шаров, помещающихся по вертикали в окно (задает размер шара относительно размера окна)

	var vx0 = 1 * a0 / T0;

	
	
	
	
    // *** Задание вычислительных параметров ***

    const fps = 50;                         // frames per second - число кадров в секунду (качеcтво отображения)
    const spf = 100;                        // steps per frame   - число шагов интегрирования между кадрами (скорость расчета)
    var dt  = 0.2 * T0 / fps; 
	var n = 40;	// шаг интегрирования 
	slider_n.value = parseInt(n);
	number_n.value = parseInt(n);
	slider_dt.value = dt;
	number_dt.value = dt;
	
	        function setN(new_n) {
                n = new_n;
	for (var z = 0; z < n; z++) {
		b[z] = [];
		a[z] = [];
	}
	for (var i = 0; i< n; i++) {
		for (var j = 0; j < n; j++) {
			c = [];
			d = [];
			c.x = w / 3 / n * i;
			c.vx = 2 * vx0 / n * j;
			d.x = w / 3 / n * i;
			d.vx = 2 * vx0 / n * j;
			b[i][j] = c;
			a[i][j] = d;
		}
	}
		context_gr.clearRect(0, 0, w * scale, h * scale);
		context_gr1.clearRect(0, 0, w * scale, h * scale);
        }
		
		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 / 3 / n * i;
			c.vx = 2 * vx0 / n * j;
			d.x = w / 3 / n * i;
			d.vx = 2 * vx0 / n * j;
			b[i][j] = c;
			a[i][j] = d;
		}
	}
		context_gr.clearRect(0, 0, w * scale, h * scale);
		context_gr1.clearRect(0, 0, w * scale, h * scale);
        }
		
	slider_n.oninput = function () {
		number_n.value = slider_n.value;
		setN(slider_n.value);
        };
		
	number_n.oninput = function () {
		slider_n.value = number_n.value;
		setN(number_n.value);
        };

	slider_dt.oninput = function () {
		number_dt.value = slider_dt.value;
		setdt(slider_dt.value);
        };
		
	number_dt.oninput = function () {
		slider_dt.value = number_dt.value;
		setdt(number_dt.value);
        };

	
	
	// Задание констант для рисования
	const scale = canvas_gr.height / Ny / a0;  // масштабный коэффициент для перехода от расчетных к экранным координатам
	
	var w = canvas_gr.width / scale;           // ширина окна в расчетных координатах
    var h = canvas_gr.height / scale;          // высота окна в расчетных координатах
	

    // -------------------------------         Выполнение программы              ------------------------------------------
	// Добавление шара
	var b = [];
	var a = [];
	for (var z = 0; z < n; z++) {
		b[z] = [];
		a[z] = [];
	}
	for (var i = 0; i< n; i++) {
		for (var j = 0; j < n; j++) {
			c = [];
			d = [];
			c.x = w / 3 / n * i;
			c.vx = 2 * vx0 / n * j;
			d.x = w / 3 / n * i;
			d.vx = 2 * vx0 / n * j;
			b[i][j] = c;
			a[i][j] = d;
		}
	}
	
	// Основной цикл программы
	setInterval(control, 1500 / fps);  // функция control вызывается с периодом, определяемым вторым параметром
	
// ---------------------------------------------------------------------------------------------------------------------
// ---------------------------------           Определение всех функций              -----------------------------------
// ---------------------------------------------------------------------------------------------------------------------
	
	// основная функция, вызываемая в программе
	function control() 
	{
        physics(); // делаем spf шагов интегрирование
		draw_gr(); // рисуем график
    }

	
    // Функция, делающая spf шагов интегрирования
    function physics() {                    // то, что происходит каждый шаг времен	
		for (var s = 1; s <= spf; s++) {
			for (var k = 0; k < n; k++) {
				for (var p = 0; p < n; p++) {	
			
			k1v = dt*(-( b[k][p].x));
			k2v = dt*(-( b[k][p].x)  + k1v/2);
			k3v = dt*(-( b[k][p].x)  + k2v/2);
			k4v = dt*(-( b[k][p].x) + 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; 
        }
				}
			}			
		}
		

	
		

	// Определение функции, рисующей график
	context_gr.fillStyle = "#3070d0"; // цвет
	context_gr.strokeStyle = "#ff0000";
	context_gr1.fillStyle = "#3070d0"; // цвет
	context_gr1.strokeStyle = "#ff0000";
	function draw_gr() 
	{
		context_gr.clearRect(0, 0, w * scale, h * scale);
		context_gr.beginPath();
		context_gr.strokeStyle = "#000000";
		
		context_gr1.clearRect(0, 0, w * scale, h * scale);
		context_gr1.beginPath();
		context_gr1.strokeStyle = "#000000";
		
		// ось
		context_gr.moveTo(w/2*scale, (h)*scale);
		context_gr.lineTo(w/2*scale, -h*scale);
		context_gr1.moveTo(w/2*scale, (h)*scale);
		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 + w/2)* 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 + w/2)* scale, (-a[l][m].vx + h / 2) * scale, 1, 0, 2 * Math.PI, false);                  //рисуем шар
				context_gr1.fill();
				context_gr1.closePath();
				
			}
		}
	
	}	
}
</syntaxhighlight>
</div>
</div>

=='''Выводы'''==
Мы наблюдаем симплектический характер метода lepfrog: метод почти сохраняет энергию. Метод Рунге-Кутты напротив не сохраняет энергию системы.

=='''Ссылки'''==
*[https://en.wikipedia.org/wiki/Leapfrog_integration Описание метода leapfrog]
*[https://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%A0%D1%83%D0%BD%D0%B3%D0%B5_%E2%80%94_%D0%9A%D1%83%D1%82%D1%82%D1%8B Описание метода Рунге-Кутты]
*[[Курсовые работы по ВМДС: 2015-2016]]
*[[Введение в механику дискретных сред]]
* [[Виртуальная лаборатория]]
<br/>
[[Category: Виртуальная лаборатория]]
[[Category: Программирование]]
[[Category: JavaScript]]