www.gusucode.com > pde 案例源码 matlab代码程序 > pde/SolvePDEWithCoefficientsInFunctionalFormExample.m
%% Solve PDE with Coefficients in Functional Form %% % This example shows how to write PDE coefficients in character form and in % functional form for 2-D geometry. %% Geometry % The geometry is a rectangle with a circular hole. Create a PDE model % container, and incorporate the geometry into the container. 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 %% PDE Coefficients % The PDE is parabolic, % % $$d \frac{\partial u}{\partial t} - \nabla \cdot \left( c \nabla u \right) + au= f$$ % % with the following coefficients: % % * d = 5 % * a = 0 % * f is a linear ramp up to 10, holds at 10, then ramps back down to 0: % % $$ % u = 10 \cdot % \left\{ % \begin{array}{ll} % 10t & 0 \leq t \leq 0.1 \\ 1 & 0.1 \leq t \leq 0.9 \\ 10-10t & 0.9 \leq t \leq 1 % \end{array} \right. % $$ % % * $c = 1 + x^2 + y^2$ % % Write a function for the |f| coefficient. % function 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 %% Boundary Conditions % The boundary conditions on the outer boundary (segments 1 through 4) are % Dirichlet, with the value $u \left( x,y \right) = t \left( x - y % \right)$, where $t$ is time. Suppose the circular boundary (segments 5 % through 8) has a generalized Neumann condition, with $q = 1$ and $g = x^2 % + y^2$. 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'); %% % The boundary conditions are the same as in % <http://www.mathworks.com/help/pde/ug/boundary-conditions-by-writing-functions.html?searchHighlight=Boundary%20Conditions%20for%20Scalar%20PDE#btj13_7 % Boundary Conditions for Scalar PDE>. That description uses the older % function form for specifying boundary conditions, which is no longer % recommended. This description uses the recommended object form. %% Initial Conditions % The initial condition is $u \left( x,y \right) = 0$ at $ t = 0$. u0 = 0; %% Mesh % Create the mesh. generateMesh(model); %% Parabolic Solution Times % Set the time steps for the parabolic solver to 50 steps from time 0 to % time 1. tlist = linspace(0,1,50); %% Solution % Solve the parabolic PDE. d = 5; a = 0; f = 'framp(t)'; c = '1 + x.^2 + y.^2'; u = parabolic(u0,tlist,model,c,a,f,d); %% % View an animation of the solution. 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 %% Alternative Coefficient Syntax % Equivalently, you can write a function for the coefficient |f| in the % syntax described % in <http://www.mathworks.com/help/pde/ug/scalar-coefficients-in-function-form.html?searchHighlight=Specify%202-D%20Scalar%20Coefficients%20in%20Function%20Form#responsive_offcanvas % Specify 2-D Scalar Coefficients in Function Form>. % % 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 %% % Call this function by setting f = @framp2; u = parabolic(u0,tlist,model,c,a,f,d); %% % You can also write a function for the coefficient |c|, though it is more % complicated than the character formulation. % function c = cfunc(p,t,u,time) % % % 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 %% % Call this function by setting c = @cfunc; u = parabolic(u0,tlist,model,c,a,f,d);