pdetoolbox-0802

Introduction to the Matlab Partial Differential Equation Toolbox------PDETOOL Introduction to the Matlab Partial Differential Equation Toolbox——PDETOOL Matlab PDE Toolbox介绍 (1) 可以求解(非)线性椭圆型PDE,以及线性抛物,双曲问题。可以求解单个方程，也可以求解方程组： (2) 自适应网格加密, 区域分裂(如Lshape区域分为3部分进行计算, matlab有演示算例) (3) 有两种方式调用偏微分方程工具箱中的函数： ​ matlab script ​ GUI(by >>pdetool) 如果采用GUI方式,一般的求解过程是： ​ Starting the PDE Toolbox 开启工具箱 ​ Specifying the Domain 指定二维求解区域 ​ Specifying the Boundary Conditions 指定边界条件 ​ Specifying the PDE 指定PDE方程(椭圆、抛物、双曲) ​ Specifying the Initial Mesh 指定初始网格 ​ Mesh Refinement 设置网格加密 ​ Solving boundary value problem and plotting solution 求解(初)边值问题并绘图 ​ Working in the MATLAB workspace (export data｛p,e,t,u｝) 对变量空间进行操作 ​ Post-Processing the Solution 后处理数值解比如tri2grid(p,t,u,x,y) ​ Visualization Commands 比如绘图命令pdemesh(p,e,t), pdesurf(p,t,u) (4) 非结构三角网格的网格数据结构 下面是一个半圆形区域的非结构网格的剖分结果： (5) 一般的椭圆问题 的系数矩阵c的输入法则： 注意c可以写成 4 -2 1是因为4 -2 -2 1是一个对称的二阶矩阵。Matlab PDE Toolbox 允许这样。类似地 4 0 0 1可写成4 1，即只给出二阶系数矩阵的对角线元素。 注意上面每个 都是一个 的矩阵, 在GUI中输入时写在一行按照 的顺序输入(因为Matlab的二维数组是按列优先存放的). (6) 非线性问题的2个例子, 非线性是指 中的函数c,a,f依赖于u，ux或者uy。 注： (7) Matlab的一些内置的区域名称以及边界条件名称(以便采用脚本形式灵活求解方程,不用GUI那样点来点去,选来选去的,这正是脚本方式的好处). ​ >> g='lshapeg'; >> [p,e,t]=initmesh(g); >> pdemesh(p,e,t) %下面左侧网格图形 ​ >> [p,e,t]=refinemesh(g,p,e,t); >> pdemesh(p,e,t) %中间图形(四分三角形均匀加密) ​ >> g='squareg'; n=10; [p,e,t]=poimesh(g,n); pdemesh(p,e,t) % poimesh产生结构网格,最右侧 一些定义好的边界条件: ​ >> g='circleg' %单位圆,圆心(0,0),半径1. >> g='squareg' %正方形,[-1,1]^2 ​ >> b=’circleb1’ %函数值0边界条件（齐次Dirichlet边界条件），单位圆区域 ​ >> b=’circleb2’ %函数值x.^2边界条件（非齐次Dirichlet边界条件） ​ >> b=’lshapeb’ %函数值0边界条件（齐次Dirichlet边界条件），L型区域 ​ >> b=’squareb4’ %三面0边界条件，最右端半正弦条件，正方形[-1,1]^2 于是对于已知真解的椭圆问题，可以测试Matlab的pdetool中线性元的精度。 以下程序段是在均匀加密网格上对u=sin(pi*x).*sin(pi*y)在[-1,1]^2 进行验证. 程序的运行结果见下图：(可见对于L2误差有2阶精度) Matlab脚本函数代码 function testrefinemesh(N) %-------- define the problem domain and boundary condition -------- g = 'squareg'; % form the square of [-1,1]^2 in x-y plane [p,e,t] = initmesh(g,'jiggle','on'); % ('lshapeg','jiggle','off') for i = 1:N [p,e,t] = refinemesh(g,p,e,t); end b = 'squareb1'; % set the zero Dirichlet boundary condition c = 1; % diffusion coefficient a = 0; f = '2*pi^2*sin(pi*x).*sin(pi*y)'; % right hand side function u = assempde(b,p,e,t,c,a,f); % global assembly pdeplot(p,e,t,'xydata',u,'zdata',u,'mesh','on')% pdemesh(p,e,t,u) %------ compute the L2 error of the numerical solution -------- x = p(1,:)'; y = p(2,:)'; uex = sin(pi*x).*sin(pi*y); % exact solution ne = max(size(t)); % total # of elements L2 = 0; for i = 1:ne xi = p(1,t(1,i)); yi = p(2,t(1,i)); xj = p(1,t(2,i)); yj = p(2,t(2,i)); xm = p(1,t(3,i)); ym = p(2,t(3,i)); m = [xi, yi, 1; xj, yj, 1; xm, ym, 1]; s = det(m)/2; L2 = L2 + ( (u(t(1,i))+u(t(2,i)) - uex(t(1,i))-uex(t(2,i)))^2/4 ... + (u(t(1,i))+u(t(3,i)) - uex(t(1,i))-uex(t(3,i)))^2/4 ... + (u(t(2,i))+u(t(3,i)) - uex(t(2,i))-uex(t(3,i)))^2/4 )*8*s/60 ... + ( ( u(t(1,i)) - uex(t(1,i)) )^2 ... + ( u(t(2,i)) - uex(t(2,i)) )^2 ... + ( u(t(3,i)) - uex(t(3,i)) )^2 )*3*s/60 ... + ( u(t(1,i))+u(t(2,i))+u(t(3,i)) - uex(t(1,i))-uex(t(2,i))-uex(t(3,i)) )^2/9 * 27*s/60; end L2_err = L2^0.5; disp(['There are : ',num2str(size(p,2)),' points in the FEM mesh']) disp(['The L2 error is : ',num2str(L2_err)]) (8) Matlab PDE Toolbox的函数名称 adaptmesh %生成自适应网格并求解PDE问题 assema %组合面积的整体贡献 assemb %组合边界的整体贡献 assempde %组合刚度矩阵与PDE问题的右端项 csgchk % checks if the solid objects in the Geometry Description matrix gd are valid csgdel % [dl1,bt1]=csgdel(dl,bt,bl) deletes the border segments in the list bl decsg % analyzes the Constructive Solid Geometry model (CSG model) that you draw dst, idst % computes the discrete sine transform of the columns of x hyperbolic %求解双曲问题 initmesh % [p,e,t]=initmesh(g) returns a triangular mesh using the geometry g jigglemesh % jiggles the mesh by adjusting the node positions. The quality of mesh normally increases parabolic %求解抛物问题 pdeadgsc % Select triangles using a relative tolerance criterion pdeadworst % Select triangles relative to the worst value pdearcl % Interpolation between parametric representation and arc length pdebound % Boundary M-file, [q,g,h,r]=pdebound(p,e,u,time) pdecgrad % The flux of a PDE solution pdecirc % Draw circle pdecont % Shorthand command for contour plot, pdecont(p,t,u) pdeeig %求解椭圆特征值问题 pdeellip % Draw ellipse, pdeellip(xc,yc,a,b,phi) pdeent % Indices of triangles neighboring a given set of triangles pdegeom % Geometry M-file pdegplot % Plot a PDE geometry, pdegplot(g) pdegrad % The gradient of a PDE solution, [ux,uy]=pdegrad(p,t,u) pdeintrp % Interpolate from node data to triangle midpoint data, ut=pdeintrp(p,t,un) pdejmps % Error estimates for adaption, pdemdlcv % Convert PDE Toolbox 1.0 Model M-files to PDE Toolbox 1.0.2 format pdemesh % Plot a PDE triangular mesh, pdemesh(p,e,t) pdenonlin % Solve nonlinear PDE problem(elliptic equation or systems) pdeplot % pdeplot(p,e,t,'PropertyName',PropertyValue,) pdepoly % Draw polygon, pdepoly(x,y) pdeprtni % Interpolate from triangle midpoint data to node data, un=pdeprtni(p,t,ut) pderect % Draw rectangle, pderect(xy) pdesdp, pdesde, pdesdt % Indices of points/edges/triangles in a set of subdomains pdesmech % Calculate structural mechanics tensor functions pdesurf % Shorthand command for surface plot, pdesurf(p,t,u) pdetool % PDE Toolbox graphical user interface (GUI) pdetrg % Triangle geometry data pdetriq % Triangle quality measure, q=pdetriq(p,t) poiasma % Boundary point matrix contributions for fast solvers of Poisson's equation poicalc % Fast solver for Poisson's equation on a rectangular grid poiindex % Indices of points in canonical ordering for rectangular grid poimesh % Make regular mesh on a rectangular geometry poisolv % Fast solution of Poisson's equation on a rectangular grid refinemesh % Refine a triangular mesh, [p1,e1,t1]=refinemesh(g,p,e,t) sptarn % Solve generalized sparse eigenvalue problem tri2grid % Interpolate from PDE triangular mesh to rectangular grid, uxy=tri2grid(p,t,u,x,y) wbound % Write boundary condition specification file, fid=wbound(bl,mn) wgeom % Write geometry specification function (8) 根据以上函数, 下面是一个脚本形式的Matlab PDE Toolbox Adaptive-FEM example: % Solve the L-shaped domain problem by using the adaptive finite % algorithm based on the greedy strategy. See "lshaped_uniform.m" % for a description of the L-shaped domain problem. alfa=0.15;beta=0.15;mexp=1; %For the a posteriori error estimate J=17; %Maximum number of iterations. %Geometry g = [2 2 2 2 2 2 0 1 1 -1 -1 0 1 1 -1 -1 0 0 0 0 1 1 -1 -1 0 1 1 -1 -1 0 1 1 1 1 1 1 0 0 0 0 0 0]; %Boundary conditions b=[1 1 1 1 1 52 48 48 49 40 120 46 94 50 43 121 ... 46 94 50 41 46 94 40 49 47 51 41 46 42 115 105 ... 110 40 50 47 51 42 97 116 97 110 50 40 121 44 120 ... 41 43 52 47 51 42 112 105 42 40 121 60 48 41 41]’; b=repmat(b,1,6); %PDE coefficients c=1; a=0; f=0; %Initial mesh [p,e,t]=initmesh(g); % Do iterative adaptive refinement, solve PDE, estimate the error. % Exact solution u ue=’(x.^2+y.^2).^(1/3).*sin(2/3*(atan2(y,x)+2*pi*(y<0)))’; % ux=-2/3*r^(-1/3)*sin(theta/3) uex=’-2/3*(x.^2+y.^2).^(-1/6).*sin(1/3*(atan2(y,x)+2*pi*(y<0)))’; % uy=2/3*r^(-1/3)*cos(theta/3) uey=’2/3*(x.^2+y.^2).^(-1/6).*cos(1/3*(atan2(y,x)+2*pi*(y<0)))’; error=[]; N_k=[]; for k=1:J+1 fprintf(’Number of triangles: %g\n’,size(t,2)) u=assempde(b,p,e,t,c,a,f); err=pdeerrH1(p,t,u,ue,uex,uey); %H^1 error error=[error err]; N_k=[N_k,size(p,2)]; %DoFs if k表

关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf