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

window.addEventListener("load", main_equations, false);
function main_equations() {

    // Предварительные установки

    var a0 = 1;                             // масштаб расстояния
    var t0 = 1;                             // масштаб времени
    var m0 = 1;                             // масштаб массы
    var T0 = 1;                             // масштаб температуры

    var k0 = 2 * Math.PI / t0;              // масштаб частоты
    var C0 = m0 * k0 * k0;                  // масштаб жесткости

    var m = 1 * m0;                 	    // масса
    var C = 1 * C0;                 	    // жесткость

    var omega = Math.sqrt(C / m);
    var c = a0 * omega;

    var tw = CONFIG.N_Wave * a0 / c / 4;
    var beta =  Math.PI * c * c * Math.sqrt(tw);

    var dx = 1 * a0;                        // шаг сетки по оси x

    var spf = CONFIG.spf;                   // steps per frame - сколько расчетов проходит за каждый кадр отображения
    var dt = CONFIG.dt * t0;                // шаг интегрирования по времени

    var c3 = init_canvas(Equations_canvas); // инициализация canvas

    // слайдер
    slider_dt.max = dt * 2;
    slider_dt.min = dt / 50;
    slider_dt.step = dt / 100;
    slider_dt.value = dt;
    slider_dt.oninput = function() {dt = parseFloat(slider_dt.value);};

    var heat, wave;
    function start_new_system() {
        heat = new Heat_model(T0, dx, dt, beta);
        wave = new Wave_model(a0, k0, m, c, dx, dt);
    }
    function start_new_system_stair() {
        var func = function(M, M_max, N, val) {return IC_stair(M, M_max, N, val,
            parseFloat(stair_x1_number.value), parseFloat(stair_x2_number.value))};
        CONFIG.IC_Heat = func;
        CONFIG.IC_Wave = func;
        start_new_system()
    }
    function start_new_system_Gauss() {
        CONFIG.IC_Heat = IC_Gauss;
        CONFIG.IC_Wave = IC_Gauss;
        start_new_system()
    }



    button_restart_stair.onclick = start_new_system_stair;
    button_restart_Gauss.onclick = start_new_system_Gauss;

    start_new_system();
    control();

    function control() {
        heat.calculate_steps(spf, dt);
        wave.calculate_steps(spf, dt);
        draw_interface(c3);
        if (checkbox_graph_01.checked) draw_graph(c3, heat.get_D(), heat.get_div(), 1, "#ff0000", "u");
        if (checkbox_graph_02.checked) draw_graph(c3, wave.get_D(), wave.get_div(), 1, "#00ff00", "u");
        requestAnimationFrame(control);
    }

    function draw_graph(c, D, D_div, mult_y, color, val) {
        // график
        c.ctx.strokeStyle = color;
        c.ctx.beginPath();
        c.ctx.moveTo(0, c.h - D[0].u * (c.h - 20));
        for (var i = 1; i < D.length - 1; i++) {
            c.ctx.lineTo(i / (D.length - 1) * c.w, c.h -  D[i][val] / D_div * mult_y * (c.h - 20));
        }
        c.ctx.stroke();
    }

    function draw_interface(c) {
        c.ctx.clearRect(0, 0, c.w, c.h);

        // вертикальная линия
        c.ctx.strokeStyle = "#bbbbbb";
        c.ctx.beginPath();
        c.ctx.moveTo(c.w / 2, 0);
        c.ctx.lineTo(c.w / 2, c.h);
        c.ctx.stroke();
    }

    function init_canvas(canvas) {
        canvas.onselectstart = function () {return false;};     // запрет выделения canvas

        var canv_obj = {};
        canv_obj.ctx = canvas.getContext("2d");                 // на context происходит рисование
        canv_obj.w = canvas.width;                              // ширина окна в расчетных координатах
        canv_obj.h = canvas.height;                             // высота окна в расчетных координатах

        return canv_obj;
    }
}


function Heat_model(T0, dx, dt, X) {
    var N = CONFIG.N_Heat + 2;                                  // количество узлов по оси x + 2 для ГУ
    var T_start = 1 * T0;                                       // начальная температура

    var T = make_mass(N);
    var IC = CONFIG.IC_Heat;
    IC(T, T_start, N, "u");                                     // начальные условия
    BC_mirrow(T, N);                                            // граничные условия

    var div = 1 * T0;                                           // делитель для функции рисования
    this.get_div = function() {return div};
    var t = dt / 2;                                             // суммарное время с начала расчета

    this.get_D = function() {return T};
    this.calculate_steps = calc_implicit_Euler_method;


    var alpha = new Array(N + 2), beta = new Array(N + 2);
    function calc_implicit_Euler_method(spf, dt) {
        for (var s = 0; s < spf; s++) {
            var A = X * dt / (dx * dx + 2 * X * dt);
            var B = dx * dx / (dx * dx + 2 * X * dt);

            alpha[N - 2] = 0;
            beta[N - 2] = T[N - 2].u;
            for (var i = N - 3; i >= 1; i--) {
                alpha[i] = A / (1 - alpha[i + 1] * A);
                beta[i] = (B * T[i + 1].u + beta[i + 1] * A) / (1 - alpha[i + 1] * A);
            }

            T[1].u = (B * T[1].u + beta[1] * A) / (1 - alpha[1] * A - A);
            for (var i = 2; i < N - 1; i++) {
                T[i].u = (A * T[i - 1].u + B * T[i].u + beta[i] * A) / (1 - alpha[i] * A);
            }

            t += dt;
        }
    }

    return this;
}
function Wave_model(a0, k0, m, c, dx, dt) {
    var N = CONFIG.N_Wave + 2;                                  // количество узлов по оси x + 2 для ГУ
    var U_start = a0 * k0 / (2 * Math.PI);

    var U = make_mass(N);
    var IC = CONFIG.IC_Wave;
    IC(U, U_start, N, "u");                                     // начальные условия
    IC_zero(U, U_start, N, "v");                                // начальные условия для скоростей
    BC_mirrow(U, N);                                            // граничные условия

    var div = 1 * a0 * k0 / (2 * Math.PI);                      // делитель для функции рисования
    this.get_div = function() {return div};
    var t = dt / 2;                                             // суммарное время с начала расчета

    this.get_D = function() {
        if (CONFIG.smooth_wave) return get_smoothed(U, N, "u", CONFIG.averaging_quality_Wave);
        else return U;
    };
    this.calculate_steps = function(spf, dt) {
        var koeff1 = c * c / (dx * dx) / m;
        for (var s = 0; s < spf; s++) {
            for (var i = 1; i < N - 1; i++) {
                U[i].v += (U[i + 1].u - 2 * U[i].u + U[i - 1].u) * koeff1 * dt;
            }
            for (var i = 1; i < N - 1; i++) {
                U[i].u += U[i].v * dt;
            }
            t += dt;
        }
    };
    return this;
}

function get_smoothed(D0, N, val, avl_quality) {
    var D = [];
    var avl = (N - 2) / avl_quality;
    var s = 1;
    var avl_half = avl / 2;
    for (var i = avl_half + 1; i < N - 1 ; i += avl) {
        var av = 0;
        for (var j = i - avl_half; j < i + avl_half; j++) {
            av += D0[j][val];
        }
        D[s] = {};
        D[s][val] = av / avl;
        s++;
    }
    // для первой точки добавляется объект, чтобы график располагался правильно по центру
    // но сама первая точка не добавляется - она не усреднена
    D[0] = {};
    return D;
}

var CONFIG = {
     fps:       60                              // количество кадров в секунду
    ,spf:       200                             // сколько раз происходит вычисление между кадрами

    ,dt: 0.005

    // начальные условия
    //      IC_zero:                все значения - нули
    //      IC_stair:               ступенька посередине экрана
    //      IC_half_stair:          половина ступеньки
    //      IC_stair_random:        ступенька, составленная из случайных значений
    //      IC_half_stair_random:   половина ступеньки из случайных значений
    //      IC_Gauss:               распределение Гаусса
//    ,IC_Heat:   IC_Gauss
//    ,IC_Wave:   IC_Gauss
    ,IC_Heat:   function(M, M_max, N, val) {return IC_stair(M, M_max, N, val, 0.33, 0.66)}
    ,IC_Wave:   function(M, M_max, N, val) {return IC_stair(M, M_max, N, val, 0.33, 0.66)}

    ,N_Heat:    1000                            // количество частиц по оси x
    ,N_Wave:    1000                            // количество частиц по оси x

    ,averaging_quality_Wave: 200                // на сколько отрезков усреднения поделить N (число должно быть делителем N)
    ,smooth_wave: true                         // сглаживать волновое уравнение
};

// Аргументы в функции IC и BC:
//      M - массив (сформированный с помощью функции make_mass данной библиотеки)
//      M_max - максимальное значение, которое должна задать функция (иногда данное значение передается,
//          но не испозьзуется, сделано для того, чтобы у всех фунций был одинаковый набор аргументов,
//          и можно было бы быстро менять используемую функцию (например, из файла настроек))
//      N - Размер массива
//      val - ключ, по которому будет организован доступ к данным. В программе можно будет получить данные через него,
//          например, так: M[0].val
//          Чтобы лучше разобраться - просмотрите код какой либо программы, использующей эту библиотеку.

/**
 * Создает специальный массив для последующего задания в нем начальных и граничных условий
 * N - Длина массива (нужно учитывать, что первый и последний элементы будут использованы для граничных условий))
 * */
function make_mass(N) {
    var M = [];
    for (var i = 1; i < N - 1; i++)
        M[i] = {};
    return M;
}

function IC_zero(M, M_max, N, val) {            // Заполняет массив нулями
    for (var i = 1; i < N - 1; i++) {
        M[i][val] = 0;
    }
}

// 0 <= start, end <= 1
// start <= end
// чтобы получить функцию вида IC_stair(M, M_max, N, val), используйте замыкание
// например: var func = function(M, M_max, N, val) {return IC_stair(M, M_max, N, val, 0.33, 0.66)}
function IC_stair(M, M_max, N, val, start, end) {   // Задает "ступеньку" высоты M_max посередине экрана
    for (var i = 1; i < N - 1; i++) {
        if (i >= N * start && i < N * end) M[i][val] = M_max;
        else M[i][val] = 0;
    }
}
function IC_half_stair(M, M_max, N, val) {      // Задает мгновенный подъем от нуля до максимального значения
    for (var i = 1; i < N - 1; i++) {
        if (i <  N / 2) M[i][val] = 0;
        else M[i][val] = M_max;
    }
}
function IC_stair_random(M, M_max, N, val) {    // Задает "ступеньку" из случайных значений от -M_max до M_max
    for (var i = 1; i < N - 1; i++) {
        if (i >= N / 3 && i < 2 * N / 3) M[i][val] = (2 * Math.random() - 1) * M_max;
        else M[i][val] = 0;
    }
}
function IC_half_stair_random(M, M_max, N, val) {// Правая половина экрана будет заполнена случайными значениями от -M_max до M_max
    for (var i = 1; i < N - 1; i++) {
        if (i <  N / 2) M[i][val] = 0;
        else M[i][val] = (2 * Math.random() - 1) * M_max;
    }
}
function IC_Gauss(M, M_max, N, val) {           // Задает распределение Гаусса.
    for (var i = 1; i < N - 1; i++) {
        var x = (i - 1) / (N - 2);
        M[i][val] = Math.exp(-Math.pow(x - 0.5, 2) * 50);
    }
}

function BC_mirrow(M, N) {                      // Зеркальные граничные условия
    M[0] = M[1];
    M[N - 1] = M[N - 2];
}
function BC_periodic(M, N) {                    // Периодические граничные условия.
    M[0] = M[N - 2];
    M[N - 1] = M[1];
}

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