Можно интерполировать решение и, в случае необходимости, его градиент на отдельных шагах, и затем построить результаты при помощи функций MATLAB®, таких как surf
mesh
quiver
, и так далее. Например, решите ту же скалярную эллиптическую задачу на L-образной мембране с нулем граничные условия Дирихле. Интерполируйте решение и его градиент, и затем постройте результаты.
Создайте модель PDE, 2D геометрию и mesh. Задайте граничные условия и коэффициенты. Решите задачу УЧП.
model = createpde; geometryFromEdges(model,@lshapeg); applyBoundaryCondition(model,'dirichlet', ... 'Edge',1:model.Geometry.NumEdges, ... 'u',0); c = 1; a = 0; f = 1; specifyCoefficients(model,'m',0,'d',0,'c',c,'a',a,'f',f); generateMesh(model,'Hmax',0.05); results = solvepde(model);
Интерполируйте решение и его градиенты к плотной сетке от-1 до 1 в каждом направлении.
v = linspace(-1,1,101); [X,Y] = meshgrid(v); querypoints = [X(:),Y(:)]'; uintrp = interpolateSolution(results,querypoints);
Постройте получившееся решение на mesh.
uintrp = reshape(uintrp,size(X)); mesh(X,Y,uintrp) xlabel('x') ylabel('y')
Интерполируйте градиенты решения сетки от-1 до 1 в каждом направлении. Постройте результат с помощью quiver
.
[gradx,grady] = evaluateGradient(results,querypoints); figure quiver(X(:),Y(:),gradx,grady) xlabel('x') ylabel('y')
Увеличьте масштаб, чтобы видеть больше деталей. Например, ограничьте область значений [-0.2,0.2]
в каждом направлении.
axis([-0.2 0.2 -0.2 0.2])
Постройте решение и градиенты на той же области значений.
figure h1 = meshc(X,Y,uintrp); set(h1,'FaceColor','g','EdgeColor','b') xlabel('x') ylabel('y') alpha(0.5) hold on Z = -0.05*ones(size(X)); gradz = zeros(size(gradx)); h2 = quiver3(X(:),Y(:),Z(:),gradx,grady,gradz); set(h2,'Color','r') axis([-0.2,0.2,-0.2,0.2])
Срез решения строит вдоль линии x = y
.
figure mesh(X,Y,uintrp) xlabel('x') ylabel('y') alpha(0.25) hold on z = linspace(0,0.15,101); Z = meshgrid(z); surf(X,X,Z') view([-20 -45 15]) colormap winter
Постройте интерполированное решение вдоль линии.
figure xq = v; yq = v; uintrp = interpolateSolution(results,xq,yq); plot3(xq,yq,uintrp) grid on xlabel('x') ylabel('y')
Интерполируйте градиенты решения вдоль той же линии и добавьте их в график решения.
[gradx,grady] = evaluateGradient(results,xq,yq);
gradx = reshape(gradx,size(xq));
grady = reshape(grady,size(yq));
hold on
quiver(xq,yq,gradx,grady)
view([-20 -45 75])