← Назад Далее →

3. Автоматическая генерация разностных схем и их паспортов

В данном разделе описывается алгоритм, примененный для автоматической генерации атласа диссипативных и дисперсионных свойств разностных схем второго порядка аппроксимации на компактных вычислительных шаблонах. Шаблон будем называть компактным, если он состоит из четырех узлов, индексы (j,k) которых удовлетворяют условиям i−1≤j≤i+1, n−1≤k≤n+1. Генерация паспорта каждой разностной схемы, состоит из следующих этапов:

  1. построение очередного компактного шаблона;
  2. построение разностной схемы, соответствующей шаблону;
  3. построение символьного выражения для разностной схемы;
  4. построение характеристического полинома для схемы;
  5. поиск интервалов устойчивости разностной схемы;
  6. нахождение корней характеристического полинома внутри интервалов устойчивости;
  7. сортировка значений найденных корней;
  8. вычисление модуля и аргумента комплексных значений корней;
  9. коррекция значений аргумента;
  10. построение графиков диссипативной и дисперсионной поверхностей.

Алгоритм реализован на языке системы символьного программирования Wolfram Mathematica 10.3. Ниже будут приведены некоторые фрагменты кода, иллюстрирующие некоторые особенности реализации алгоритма. Комментарии в коде расположены после операторов, к которым они относятся.

Основной цикл программы выглядит следующим образом:

initial = Table[If[i < = 4, 1, 0], {i, 9}]; (* начальная последовательность *)
templates = Permutations@initial; (* список всех возможных перестановок *)
For[iTask = 1, iTask < Length[templates], iTask++,    
    task = templates; (* очередной шаблон *)
    sсh = BuildSсheme[task]; (*построение схемы по шаблону *)
    sсhout = BuildShOut[sсh]; (* символьное представление схемы *)
    poly = BuildPoly[sсh]; (* характеристический полином *)
    limits = BuildLimits[poly]; (* интервалы устойчивости *)
    For[ilim = 1, ilim < Length[ilim], ilim++, (* цикл по интервалам устойчивости *)
        limit = limits;   
        rmin = limit; (* нижний предел интервала *)
        rmax = limit; (* верхний предел интервала *)
        If[rmin == rmax, Continue[]]; (* пропуск вырожденного интервала *)
        roots = BuildRoots[sh, rmin, rmax]; (* нахождение и сортировка корней *)
        sol = BuildDis2[roots, rmin, rmax]; (* вычисление модуля и аргумента *)
        dis2 = CheckSolution[sol, rmin, rmax]; (*коррекция значений *)
        surfs = BuildGraphs[dis2, rmin, rmax]; (*построение графиков *)
        OutGraphs[surfs, rmin, rmax]; (* вывод в файл *)
    ];
];

← Назад Далее →