Численное решение уравнения теплопроводности и волнового уравнения — различия между версиями

  1 window.addEventListener("load", main_equations, false);
  2 function main_equations() {
  3 
  4     // Предварительные установки
  5 
  6     var a0 = 1;                             // масштаб расстояния
  7     var t0 = 1;                             // масштаб времени
  8     var m0 = 1;                             // масштаб массы
  9     var T0 = 1;                             // масштаб температуры
 10 
 11     var k0 = 2 * Math.PI / t0;              // масштаб частоты
 12     var C0 = m0 * k0 * k0;                  // масштаб жесткости
 13 
 14     var m = 1 * m0;                 	    // масса
 15     var C = 1 * C0;                 	    // жесткость
 16 
 17     var omega = Math.sqrt(C / m);
 18     var c = a0 * omega;
 19 
 20     var tw = CONFIG.N_Wave * a0 / c / 4;
 21     var beta =  Math.PI * c * c * Math.sqrt(tw);
 22 
 23     var dx = 1 * a0;                        // шаг сетки по оси x
 24 
 25     var spf = CONFIG.spf;                   // steps per frame - сколько расчетов проходит за каждый кадр отображения
 26     var dt = CONFIG.dt * t0;                // шаг интегрирования по времени
 27 
 28     var c3 = init_canvas(Equations_canvas); // инициализация canvas
 29 
 30     // слайдер
 31     slider_dt.max = dt * 2;
 32     slider_dt.min = dt / 50;
 33     slider_dt.step = dt / 100;
 34     slider_dt.value = dt;
 35     slider_dt.oninput = function() {dt = parseFloat(slider_dt.value);};
 36 
 37     var heat, wave;
 38     function start_new_system() {
 39         heat = new Heat_model(T0, dx, dt, beta);
 40         wave = new Wave_model(a0, k0, m, c, dx, dt);
 41     }
 42     function start_new_system_stair() {
 43         var func = function(M, M_max, N, val) {return IC_stair(M, M_max, N, val,
 44             parseFloat(stair_x1_number.value), parseFloat(stair_x2_number.value))};
 45         CONFIG.IC_Heat = func;
 46         CONFIG.IC_Wave = func;
 47         start_new_system()
 48     }
 49     function start_new_system_Gauss() {
 50         CONFIG.IC_Heat = IC_Gauss;
 51         CONFIG.IC_Wave = IC_Gauss;
 52         start_new_system()
 53     }
 54 
 55 
 56 
 57     button_restart_stair.onclick = start_new_system_stair;
 58     button_restart_Gauss.onclick = start_new_system_Gauss;
 59 
 60     start_new_system();
 61     control();
 62 
 63     function control() {
 64         heat.calculate_steps(spf, dt);
 65         wave.calculate_steps(spf, dt);
 66         draw_interface(c3);
 67         if (checkbox_graph_01.checked) draw_graph(c3, heat.get_D(), heat.get_div(), 1, "#ff0000", "u");
 68         if (checkbox_graph_02.checked) draw_graph(c3, wave.get_D(), wave.get_div(), 1, "#00ff00", "u");
 69         requestAnimationFrame(control);
 70     }
 71 
 72     function draw_graph(c, D, D_div, mult_y, color, val) {
 73         // график
 74         c.ctx.strokeStyle = color;
 75         c.ctx.beginPath();
 76         c.ctx.moveTo(0, c.h - D[0].u * (c.h - 20));
 77         for (var i = 1; i < D.length - 1; i++) {
 78             c.ctx.lineTo(i / (D.length - 1) * c.w, c.h -  D[i][val] / D_div * mult_y * (c.h - 20));
 79         }
 80         c.ctx.stroke();
 81     }
 82 
 83     function draw_interface(c) {
 84         c.ctx.clearRect(0, 0, c.w, c.h);
 85 
 86         // вертикальная линия
 87         c.ctx.strokeStyle = "#bbbbbb";
 88         c.ctx.beginPath();
 89         c.ctx.moveTo(c.w / 2, 0);
 90         c.ctx.lineTo(c.w / 2, c.h);
 91         c.ctx.stroke();
 92     }
 93 
 94     function init_canvas(canvas) {
 95         canvas.onselectstart = function () {return false;};     // запрет выделения canvas
 96 
 97         var canv_obj = {};
 98         canv_obj.ctx = canvas.getContext("2d");                 // на context происходит рисование
 99         canv_obj.w = canvas.width;                              // ширина окна в расчетных координатах
100         canv_obj.h = canvas.height;                             // высота окна в расчетных координатах
101 
102         return canv_obj;
103     }
104 }
105 
106 
107 function Heat_model(T0, dx, dt, X) {
108     var N = CONFIG.N_Heat + 2;                                  // количество узлов по оси x + 2 для ГУ
109     var T_start = 1 * T0;                                       // начальная температура
110 
111     var T = make_mass(N);
112     var IC = CONFIG.IC_Heat;
113     IC(T, T_start, N, "u");                                     // начальные условия
114     BC_mirrow(T, N);                                            // граничные условия
115 
116     var div = 1 * T0;                                           // делитель для функции рисования
117     this.get_div = function() {return div};
118     var t = dt / 2;                                             // суммарное время с начала расчета
119 
120     this.get_D = function() {return T};
121     this.calculate_steps = calc_implicit_Euler_method;
122 
123 
124     var alpha = new Array(N + 2), beta = new Array(N + 2);
125     function calc_implicit_Euler_method(spf, dt) {
126         for (var s = 0; s < spf; s++) {
127             var A = X * dt / (dx * dx + 2 * X * dt);
128             var B = dx * dx / (dx * dx + 2 * X * dt);
129 
130             alpha[N - 2] = 0;
131             beta[N - 2] = T[N - 2].u;
132             for (var i = N - 3; i >= 1; i--) {
133                 alpha[i] = A / (1 - alpha[i + 1] * A);
134                 beta[i] = (B * T[i + 1].u + beta[i + 1] * A) / (1 - alpha[i + 1] * A);
135             }
136 
137             T[1].u = (B * T[1].u + beta[1] * A) / (1 - alpha[1] * A - A);
138             for (var i = 2; i < N - 1; i++) {
139                 T[i].u = (A * T[i - 1].u + B * T[i].u + beta[i] * A) / (1 - alpha[i] * A);
140             }
141 
142             t += dt;
143         }
144     }
145 
146     return this;
147 }
148 function Wave_model(a0, k0, m, c, dx, dt) {
149     var N = CONFIG.N_Wave + 2;                                  // количество узлов по оси x + 2 для ГУ
150     var U_start = a0 * k0 / (2 * Math.PI);
151 
152     var U = make_mass(N);
153     var IC = CONFIG.IC_Wave;
154     IC(U, U_start, N, "u");                                     // начальные условия
155     IC_zero(U, U_start, N, "v");                                // начальные условия для скоростей
156     BC_mirrow(U, N);                                            // граничные условия
157 
158     var div = 1 * a0 * k0 / (2 * Math.PI);                      // делитель для функции рисования
159     this.get_div = function() {return div};
160     var t = dt / 2;                                             // суммарное время с начала расчета
161 
162     this.get_D = function() {
163         if (CONFIG.smooth_wave) return get_smoothed(U, N, "u", CONFIG.averaging_quality_Wave);
164         else return U;
165     };
166     this.calculate_steps = function(spf, dt) {
167         var koeff1 = c * c / (dx * dx) / m;
168         for (var s = 0; s < spf; s++) {
169             for (var i = 1; i < N - 1; i++) {
170                 U[i].v += (U[i + 1].u - 2 * U[i].u + U[i - 1].u) * koeff1 * dt;
171             }
172             for (var i = 1; i < N - 1; i++) {
173                 U[i].u += U[i].v * dt;
174             }
175             t += dt;
176         }
177     };
178     return this;
179 }
180 
181 function get_smoothed(D0, N, val, avl_quality) {
182     var D = [];
183     var avl = (N - 2) / avl_quality;
184     var s = 1;
185     var avl_half = avl / 2;
186     for (var i = avl_half + 1; i < N - 1 ; i += avl) {
187         var av = 0;
188         for (var j = i - avl_half; j < i + avl_half; j++) {
189             av += D0[j][val];
190         }
191         D[s] = {};
192         D[s][val] = av / avl;
193         s++;
194     }
195     // для первой точки добавляется объект, чтобы график располагался правильно по центру
196     // но сама первая точка не добавляется - она не усреднена
197     D[0] = {};
198     return D;
199 }

 1 var CONFIG = {
 2      fps:       60                              // количество кадров в секунду
 3     ,spf:       200                             // сколько раз происходит вычисление между кадрами
 4 
 5     ,dt: 0.005
 6 
 7     // начальные условия
 8     //      IC_zero:                все значения - нули
 9     //      IC_stair:               ступенька посередине экрана
10     //      IC_half_stair:          половина ступеньки
11     //      IC_stair_random:        ступенька, составленная из случайных значений
12     //      IC_half_stair_random:   половина ступеньки из случайных значений
13     //      IC_Gauss:               распределение Гаусса
14 //    ,IC_Heat:   IC_Gauss
15 //    ,IC_Wave:   IC_Gauss
16     ,IC_Heat:   function(M, M_max, N, val) {return IC_stair(M, M_max, N, val, 0.33, 0.66)}
17     ,IC_Wave:   function(M, M_max, N, val) {return IC_stair(M, M_max, N, val, 0.33, 0.66)}
18 
19     ,N_Heat:    1000                            // количество частиц по оси x
20     ,N_Wave:    1000                            // количество частиц по оси x
21 
22     ,averaging_quality_Wave: 200                // на сколько отрезков усреднения поделить N (число должно быть делителем N)
23     ,smooth_wave: true                         // сглаживать волновое уравнение
24 };

 1 // Аргументы в функции IC и BC:
 2 //      M - массив (сформированный с помощью функции make_mass данной библиотеки)
 3 //      M_max - максимальное значение, которое должна задать функция (иногда данное значение передается,
 4 //          но не испозьзуется, сделано для того, чтобы у всех фунций был одинаковый набор аргументов,
 5 //          и можно было бы быстро менять используемую функцию (например, из файла настроек))
 6 //      N - Размер массива
 7 //      val - ключ, по которому будет организован доступ к данным. В программе можно будет получить данные через него,
 8 //          например, так: M[0].val
 9 //          Чтобы лучше разобраться - просмотрите код какой либо программы, использующей эту библиотеку.
10 
11 /**
12  * Создает специальный массив для последующего задания в нем начальных и граничных условий
13  * N - Длина массива (нужно учитывать, что первый и последний элементы будут использованы для граничных условий))
14  * */
15 function make_mass(N) {
16     var M = [];
17     for (var i = 1; i < N - 1; i++)
18         M[i] = {};
19     return M;
20 }
21 
22 function IC_zero(M, M_max, N, val) {            // Заполняет массив нулями
23     for (var i = 1; i < N - 1; i++) {
24         M[i][val] = 0;
25     }
26 }
27 
28 // 0 <= start, end <= 1
29 // start <= end
30 // чтобы получить функцию вида IC_stair(M, M_max, N, val), используйте замыкание
31 // например: var func = function(M, M_max, N, val) {return IC_stair(M, M_max, N, val, 0.33, 0.66)}
32 function IC_stair(M, M_max, N, val, start, end) {   // Задает "ступеньку" высоты M_max посередине экрана
33     for (var i = 1; i < N - 1; i++) {
34         if (i >= N * start && i < N * end) M[i][val] = M_max;
35         else M[i][val] = 0;
36     }
37 }
38 function IC_half_stair(M, M_max, N, val) {      // Задает мгновенный подъем от нуля до максимального значения
39     for (var i = 1; i < N - 1; i++) {
40         if (i <  N / 2) M[i][val] = 0;
41         else M[i][val] = M_max;
42     }
43 }
44 function IC_stair_random(M, M_max, N, val) {    // Задает "ступеньку" из случайных значений от -M_max до M_max
45     for (var i = 1; i < N - 1; i++) {
46         if (i >= N / 3 && i < 2 * N / 3) M[i][val] = (2 * Math.random() - 1) * M_max;
47         else M[i][val] = 0;
48     }
49 }
50 function IC_half_stair_random(M, M_max, N, val) {// Правая половина экрана будет заполнена случайными значениями от -M_max до M_max
51     for (var i = 1; i < N - 1; i++) {
52         if (i <  N / 2) M[i][val] = 0;
53         else M[i][val] = (2 * Math.random() - 1) * M_max;
54     }
55 }
56 function IC_Gauss(M, M_max, N, val) {           // Задает распределение Гаусса.
57     for (var i = 1; i < N - 1; i++) {
58         var x = (i - 1) / (N - 2);
59         M[i][val] = Math.exp(-Math.pow(x - 0.5, 2) * 50);
60     }
61 }
62 
63 function BC_mirrow(M, N) {                      // Зеркальные граничные условия
64     M[0] = M[1];
65     M[N - 1] = M[N - 2];
66 }
67 function BC_periodic(M, N) {                    // Периодические граничные условия.
68     M[0] = M[N - 2];
69     M[N - 1] = M[1];
70 }

 1 <!DOCTYPE html>
 2 <html>
 3 <head>
 4     <meta charset="UTF-8" />
 5     <title>Equations</title>
 6     <script src="IC_BC_LIB_v1.js"></script>
 7     <script src="Config.js"></script>
 8     <script src="Equations.js"></script>
 9 </head>
10 <body>
11     <canvas id="Equations_canvas" width="1000" height="500" style="border:1px solid #000000;"></canvas><br>
12     Скорость расчета: <input type="range" id="slider_dt" style="width: 150px;"><br><br>
13     <input type="checkbox" id="checkbox_graph_01" checked/><font color="#ff0000" size="5"><B>—</B></font> Уравнение теплопроводности<br>
14     <input type="checkbox" id="checkbox_graph_02" checked/><font color="#00ff00" size="5"><B>—</B></font> Волновое уравнение<br><br>
15     <input type="button" id="button_restart_Gauss" value="Рестарт Гаусс"/><br><br>
16     <input type="button" id="button_restart_stair" value="Рестарт ступенька"/> (0 <= <b>x1</b>, <b>x2</b> <= 1; <b>x1</b> < <b>x2</b>)<br>
17     x1 = <input type="number" id="stair_x1_number" min=0 max=1 step=0.001 value="0.33"/><br>
18     x2 = <input type="number" id="stair_x2_number" min=0 max=1 step=0.001 value="0.66"/>
19 </body>
20 </html>