使用泊松方程测试MMS(Method of Manufactured Solutions)方法适用性

%===================================================% TEST MMS (Method of Manufactured Solutions)%===================================================functionpoissonMMS()%======================================% u=cn*(x-x0)-sn*(y-y0)% v=sn*(x-x0)+cn*(y-y0)% T(x,y)=exp(-au^2-bv^2)%======================================theta=pi/4;p.a=10;p.b=100;p.cn=cos(theta);p.sn=sin(theta);p.x0=0.5;p.y0=0.5;model=createpde();%====================================% Geometry: [0,1]x[0,1]%====================================gd=[3;4;0;1;1;0;0;0;1;1];% [3; nsides; x1,x2,x3,x4; y1,y2,y3,y4]sf='R1';ns=char('R1')';dl=decsg(gd,sf,ns);geometryFromEdges(model,dl);figure;pdegplot(model,'EdgeLabels','on');axis equal;%=======================% build mesh%=======================generateMesh(model,'Hmax',0.025);figure;pdemesh(model);%============================================================% PDE Coefficents: -∇²u = f → m=0, d=0, c=1, a=0, f = rhs%============================================================specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,...'f',@(location,state)rhs(location,p));%==============================================================% Dirichlet BCs%===============================================================applyBoundaryCondition(model,'dirichlet','Edge',1:4,...'u',@(location,state)Tfun(location,p));%================================================================% Solve%================================================================results=solvepde(model);u=results.NodalSolution;%=================================================================%plot%=================================================================%pdeplot(model, 'XYData', u, 'Contour', 'on');%title('Solution of Poisson Equation');figure;pdeplot(model,'XYData',u,'ZData',u);title('Steady-State Temperature Distribution');xlabel('x');ylabel('y');colorbar;end%=====================================% T(x,y)=exp(-au^2-bv^2)%======================================functionT=Tfun(location,p)cn=p.cn;sn=p.sn;x0=p.x0;y0=p.y0;a=p.a;b=p.b;x=location.x;y=location.y;u=cn*(x-x0)-sn*(y-y0);v=sn*(x-x0)+cn*(y-y0);T=exp(-a*u.^2-b*v.^2);end%================================================% -△T=Residual(x,y)=(2(a+b)-4(au)^2-4(bv)^2)*T%===============================================functionres=rhs(location,p)x=location.x;y=location.y;cn=p.cn;sn=p.sn;x0=p.x0;y0=p.y0;a=p.a;b=p.b;u=cn*(x-x0)-sn*(y-y0);v=sn*(x-x0)+cn*(y-y0);f=2*(a+b)-4*((a*u).^2+(b*v).^2);res=f.*exp(-a*u.^2-b*v.^2);end