ЭТОТ ПЭЙДЖ ОПИСЫВАЕТ УСТАРЕВШИЙ РАБОЧИЙ ПРОЦЕСС. Новые возможности не могут быть совместимы с устаревшим рабочим процессом. Для соответствующего шага в рекомендуемом рабочем процессе смотрите рекомендуемые примеры на странице PDE Coefficients.
Этот пример показывает, как написать коэффициенты УЧП в символьной форме и в функциональной форме для 2D геометрии.
Геометрия является прямоугольником с круговой дырой. Создайте контейнер модели PDE и включите геометрию в контейнер.
model = createpde(1); % Rectangle is code 3, 4 sides, % followed by x-coordinates and then y-coordinates R1 = [3,4,-1,1,1,-1,-.4,-.4,.4,.4]'; % Circle is code 1, center (.5,0), radius .2 C1 = [1,.5,0,.2]'; % Pad C1 with zeros to enable concatenation with R1 C1 = [C1;zeros(length(R1)-length(C1),1)]; geom = [R1,C1]; % Names for the two geometric objects ns = (char('R1','C1'))'; % Set formula sf = 'R1 - C1'; % Create geometry gd = decsg(geom,sf,ns); % Include the geometry in the model geometryFromEdges(model,gd); % View geometry pdegplot(model,'EdgeLabels','on') xlim([-1.1 1.1]) axis equal
УЧП является параболическим,
со следующими коэффициентами:
d = 5
a = 0
f является линейным подъемом до 10, содержит в 10, затем сползает, отступают к 0:
Запишите функцию для коэффициента f. функционируйте f = framp (t)
if t <= 0.1
f = 10*t;
elseif t <= 0.9
f = 1;
else
f = 10-10*t;
end
f = 10*f;
endГраничными условиями на внешней границе (сегменты 1 - 4) является Дирихле со значением
, где
время. Предположим, что круговой контур (сегменты 5 - 8) имеет обобщенное Нейманово условие, с
и
.
myufun = @(region,state)state.time*(region.x - region.y); mygfun = @(region,state)(region.x.^2 + region.y.^2); applyBoundaryCondition(model,'edge',1:4,'u',myufun,'Vectorized','on'); applyBoundaryCondition(model,'edge',5:8,'q',1,'g',mygfun,'Vectorized','on');
Граничные условия эквивалентны в Граничных условиях для Скалярного УЧП. То описание использует более старую функциональную форму для определения граничных условий, который больше не рекомендуется. Это описание использует рекомендуемую объектную форму.
Начальное условие
на уровне $ t = 0$.
u0 = 0;
Создайте mesh.
generateMesh(model,'GeometricOrder','linear');
Установите временные шаги для параболического решателя к 50 шагам со времени 0 ко времени 1.
tlist = linspace(0,1,50);
Решите параболический УЧП.
d = 5; a = 0; f = 'framp(t)'; c = '1 + x.^2 + y.^2'; u = parabolic(u0,tlist,model,c,a,f,d);
120 successful steps 9 failed attempts 260 function evaluations 1 partial derivatives 35 LU decompositions 259 solutions of linear systems
Просмотрите анимацию решения.
for tt = 1:size(u,2) % number of steps pdeplot(model,'XYData',u(:,tt),'ZData',u(:,tt),'ColorMap','jet') axis([-1 1 -1/2 1/2 -1.5 1.5 -1.5 1.5]) % use fixed axis title(['Step ' num2str(tt)]) view(-45,22) drawnow pause(.1) end
Эквивалентно, можно записать функцию для коэффициента, f в синтаксисе, описанном в, Задает 2D Скалярные Коэффициенты в Функциональной Форме.
function f = framp2(p,t,u,time)
if time <= 0.1
f = 10*time;
elseif time <= 0.9
f = 1;
else
f = 10 - 10*time;
end
f = 10*f;
endВызовите эту функцию путем установки
f = @framp2; u = parabolic(u0,tlist,model,c,a,f,d);
120 successful steps 9 failed attempts 260 function evaluations 1 partial derivatives 35 LU decompositions 259 solutions of linear systems
Можно также записать функцию для коэффициента c, хотя это более сложно, чем символьная формулировка. функционируйте c = cfunc (p, t, u, время)
% Triangle point indices it1 = t(1,:); it2 = t(2,:); it3 = t(3,:);
% Find centroids of triangles xpts = (p(1,it1) + p(1,it2) + p(1,it3))/3; ypts = (p(2,it1) + p(2,it2) + p(2,it3))/3;
c = 1 + xpts.^2 + ypts.^2; end
Вызовите эту функцию путем установки
c = @cfunc; u = parabolic(u0,tlist,model,c,a,f,d);
120 successful steps 9 failed attempts 260 function evaluations 1 partial derivatives 35 LU decompositions 259 solutions of linear systems