Сравнение различных способов расчета градиента скалярной функции на регулярных, неортогональных сетках

Получим разбиение области Ω на треугольные элементы — триангуляция области Ω. Триангуляция такого типа называется триангуляцией Фридрихса-Келлера.Будем искать приближенное решение данного уравнения как функцию , равную нулю на границе (краевое условие), непрерывную на области Ω и линейную на каждом полученном элементе триангуляции.Функцию можно представить в следующем виде:Подставив функцию в первое уравнение, осуществив преобразования и вынос констант из под знака интеграла, сведем задачу для каждой базисной функции к подсчету интегралов вида:Значение интеграла может быть не нулевым лишь в том случае, если базисные функции под знаком интеграла имеют непустую общую область определения. По построению, каждый элемент имеет три вершины. Вершина может быть общей максимально для 6 треугольников.После подсчетов интеграла уравнение с номером будет выглядеть следующим образом:гдеи при достаточно малом h:3.3. РеализациянаЭВМfunction[x,y,U]=puass_2d_mke(c,a)% CoordinateslowerLeft = [0 ,-1];lowerRight = [1 , -1];upperRight = [1 , 1];upperLeft = [0 , 1];% Geometry matrixS = [3,4 lowerLeft(1), lowerRight(1), upperRight(1), upperLeft(1), ... lowerLeft(2), lowerRight(2), upperRight(2), upperLeft(2)]; gdm = S';% Namesns = 'S';% Set formula sf = 'S';% Invoke decsgg = decsg(gdm,ns,sf');% Import g into model using geometryFromEdges.model = createpde;f1 = 'sin(y.^2)';f2 = 'cos(3.*y)';f3 = '10*sin(x.^2)';f4 = '10*sin(6.*x)';applyBoundaryCondition(model,'dirichlet','Edge',1,'r',@f1);applyBoundaryCondition(model,'dirichlet','Edge',1,'r',@f2);applyBoundaryCondition(model,'neumann','Face',2,'g',@f3,'q',0);applyBoundaryCondition(model,'neumann','Face',2,'g',@f4,'q',0);geometryFromEdges(model,g);f = 'exp(x.*(-1))+exp(y.*(-1))';[p,e,t]= initmesh(g,'Hmax',0.1); % set the meshbc='circleb1' ; % set the boundary conditionsu = assempde(bc,p,e,t,c,a,f);figure ; % plot the meshpdemesh(p,e,t),axis equalxlabel('x'), ylabel('y')title([ 'Number of mesh points np=' int2str(size(p,2)) ' ,triangles nt=' int2str(size(t,2))])figure; % plot the error .pdesurf( p , t , u)colormap ('hsv')xlabel('y')ylabel('x')zlabel('U')function bcMatrix = f1(location,state)bcMatrix = sin(location.y.^2)function bcMatrix = f2(location,state)bcMatrix = cos(3*location.y)function bcMatrix = f3(location,state)bcMatrix = 10*sin(location.x.^2)function bcMatrix = f4(location,state)bcMatrix = 10*sin(6*location.x)Рисунок 4. Результат работы скрипта.5. Сравнение методовПоскольку, градиент по сути, является прямым следствием решения уравнения пуассона, ввиду следующего уравнения:То будет справедливым, что зависимость погрешности расчёта каждого метода, для получения частных решений как это было сделано выше, эквивалентна погрешностям каждого из методов.Рисунок 5. Графики сходимости для методов.Как видим, метод конечных элементов сходится быстрее чем, метод конечных разностей, и является более точным. Поэтому в случае малых требований по точности, при достаточно низком шаге сетки, лучше использовать метод конечных разностей. поскольку при его использовании ЭВМ будет быстрее и проще считать. Наоборот, в случае повышенной точности, следует использовать метод конечных элементов.Погрешность обоих методов возрастает с при увеличении шага по близкой к экспоненциальной зависимости. Однако стоит отметить, что у метода конечных элементов, на порядок медленнее падает погрешность при увеличении шага сетки, в виду того, что сам принцип построения сетки другой(по два триангулярных пространства на одну клетку). Точное решение, по сути является абсолютно гладкой поверхностью с формой полученной по результатам работы, что очевидно показывает что второй метод эффективнее, ввиду его большей близости к гладкому пространству, из-за триангулярных сеток.ЗаключениеВ ходе работы, были исследованы такие численные методы расчёта дифф. уравнений, как метод конечных элементов и метод конечных разностей, на примере решения дифференциального уравнения Пуассона для расчёта градиента этого уравнения на регулярных, неортоганальных сетках.Расчёты проводились на ЭВМ, путём использования пакета мат. моделирования MatLab, в ходе чего были сделаны заключения касаемо данных методов.В случае использования МКЭ, выделены следующие плюсы: метод является проекционным(то есть устойчив), решение сразу представляет собой функцию и значения в любой точке могут быть вычислены сразу.В случае использования МКР:для простых задач построение разностной схемы выполняется быстрее.