偏微分方程学习

1、椭圆形方程差分方法应用带微分余项的Taylor公式可得6665544'''3''2'7201202462hOcuhcuhcuhcuhcuhchucuhcu6665544'''3''2'7201202462hOcuhcuhcuhcuhcuhchucuhcu将上两式相加，得66644''23601221hOcuhcuhcuhcucuhcuh使用标准的二阶中心差分格式，可得2,1,,1,2xuuuujijijijix6664436012xOuxuxuxxxx2,1,,1,2yuuuujijijijiy6664436012yOuyuyuyyyy我们可以分别应用符号四阶紧致逼近算子到方程1中的二阶导xxu和yyu。离散的二维泊松方程可以表示为如下形式6424121222122121121Oyxfuyuxyyxx1和2用来表示在理查德森外推过程中忽略掉的误差，6O表示六阶精度66yxO。应用符号算子，设等于0，略去高阶项，方程可以写成42222422222222221211121121121121OfyxOfyxuxuyyxyxyxxy如果我们设置网格比为yx/，则方程可以写成1,1,,1,2,,1,11,1,1,11,11,11,182jijijijijijijijijijijijijiffffxauuubuucuuuud系数是2110a，25b，152c，2/12d.当网格比为1时，差分格式为如下形式1,1,,1,2,,1,11,1,1,11,11,11,1822044jijijijijijijijijijijijijiffffxuuuuuuuuu用矩阵可表示为1-4-1-4-204-1-4-1-由于该差分格式是九点紧致格式，四阶收敛，所以误差可以表示为)()(4hOhUAeh)16()2(42hOhUAeh这里he表示在步长为h的情况下的误差，A为精确解，)(hU表示在步长为h的情况下的数值解。由此可以看出，当步长降低一倍时，误差比值为16，这个可以有效的检验计算结果的准确性。数值算例,10,10ysin12yxeux,sin,0yyu,sin,1yeyu,10y,00,xu,01,xu.10x该问题的精确解为yeuxsin.表1给出了不同步长取值下的误差，误差比值及计算时间。表1计算误差表h最大误差误差比值运行时间0.14.8242e-5－0.0665280.053.0074e-616.04110.0457960.0251.8818e-715.98160.0534910.01251.1759e-816.00240.154051误差比值接近16，因此可以看出计算结果是比较准确的。图1、图2、图3分别为步长0.05时的精确解曲面，数值解曲面和误差曲面。图1精确解曲面图2步长为0.05时的数值解曲面图3步长为0.05时的误差曲面MATLAB程序function[error]=jiudian(h)%UNTITLEDSummaryofthisfunctiongoeshere%Detailedexplanationgoesheretic;n=1/h;fori=1:n+1%boundaryconditionu(1,i)=sin(pi*h*(i-1));u(n+1,i)=exp(1)*sin(pi*h*(i-1));u(i,1)=0;u(i,n+1)=0;end%coefficientmatrixAB=sparse(20*eye(n-1));D=sparse(-4*eye(n-1));fori=1:n-2B(i,i+1)=-4;B(i+1,i)=-4;endfori=1:n-2D(i,i+1)=-1;D(i+1,i)=-1;enda=cell(n-1);fori=1:n-1forj=1:n-1a{i,j}=sparse(zeros(n-1));endendfori=1:n-1a(i,i)={B};endfori=1:n-2a{i,i+1}=D;a{i+1,i}=D;end%CellArraytoMatrixA=cell2mat(a);A=sparse(A);%rightvectorKfori=1:n+1forj=1:n+1f(i,j)=(pi*pi-1)*exp((i-1)*h)*sin(pi*(j-1)*h);endendfori=2:nforj=2:ng(i-1,j-1)=8*f(i,j)+f(i,j+1)+f(i,j-1)+f(i+1,j)+f(i-1,j);endendk=zeros(n-1);fori=1:n-1k(1,i)=4*u(1,i+1)+u(1,i)+u(1,i+2)+k(1,i);k(i,1)=4*u(i+1,1)+u(i,1)+u(i+2,1)+k(i,1);k(n-1,i)=4*u(n+1,i+1)+u(n+1,i)+u(n+1,i+2)+k(n-1,i);k(i,n-1)=4*u(i+1,n+1)+u(i,n+1)+u(i+2,n+1)+k(i,n-1);endK=h*h/2*g+k;K=reshape(K,(n-1)*(n-1),1);%solveXX=A\K;fori=2:nforj=2:nu(j,i)=X((n-1)*(i-2)+(j-1));endend%plotX,ex=0:h:1;y=0:h:1;subplot(2,2,1);surf(x,y,u);xlabel('x');ylabel('y');zlabel('u');title('差分解')fori=1:n+1forj=1:n+1p(i,j)=exp((i-1)*h)*sin(pi*(j-1)*h);endendsubplot(2,2,2);surf(x,y,p);xlabel('x');ylabel('y');zlabel('p');title('精确解');fori=1:n+1forj=1:n+1e(i,j)=abs(u(i,j)-p(i,j));endendsubplot(2,2,3);surf(x,y,e);xlabel('x');ylabel('y');zlabel('e');title('误差');error=norm(e,inf);toc;end2、有限元方法解椭圆形问题用Ritz-Galerkin方法求边值问题,21sin21cos2222222rrruyx,，,0u,,yx其中1|,22yxyx，222yxr。该边值问题的精确解为21sin,2ryxu。由于关于原点是对称的，方程右端函数仅依赖于2r。在10H上取一簇奇函数,,12jjrrwyxnj,...,2,1.为使01rj，取,12rrw因而,1,122jjrryxnj,...,2,1.以n,...,,21为基底张成n维子空间nV，nV中任一元素可表示为.1,,11221njjjnjjjnrcryxcyxu计算得,1,'yxrryxii,,,''''rryxyxjiji,2,,,10drdrddrdrdxdyyxyxajijiji,,1nji102222122,21sin21cos212,,,drrrrrrdxdyyxfyxfriiini1，将上述两式代入到Ritz-Galerkin方程组，计算即可。表2给出了取不同n值下的误差。图4为精确解曲面，图5和图6分别为n=4时的数值解和误差曲面，图7和图8分别为n=8时的数值解和误差曲面。表2不同n值的误差n误差n误差10.273284.5916e-0920.043791.9009e-1030.0050107.1185e-1244.3484e-04112.4336e-1353.1182e-05121.0436e-1461.8899e-06134.4409e-1679.9195e-08144.4409e-16图4n=4精确解曲面图5n=4的数值解曲面图6n=4的误差曲面图7N=8的数值解曲面图8N=8的误差曲面MATLAB程序function[error]=yxyty(n)symsr;fori=1:np(i)=(1-r^2)*r^(2*(i-1));d(i)=diff(p(i));endf=2*pi*cos(pi*(1-r^2)/2)+pi*pi*r*r*sin(pi*(1-r^2)/2);fori=1:nforj=1:na(i,j)=2*pi*int(r*d(i)*d(j),0,1);endendfori=1:nb(i)=2*pi*int(r*p(i)*f,0,1);endb=b';c=a\b;c=double(c);u1=0;fori=1:nu1=c(i)*r^(2*i-2)+u1;endun=(1-r^2)*u1;%xa=-1:0.05:1;%ya=-1:0.05:1;%[x,y]=meshgrid(xa,ya);%ra=(xa.^2+ya.^2).^0.5;%zz=subs(un,r,ra);%surf(x,y,zz);h=0.1;fori=-1:h:1forj=-1:h:1ra=(i^2+j^2)^0.5;a=(i+1)/h+1;a=int8(a);b=(j+1)/h+1;b=int8(b);ifra1z(a,b)=subs(un,r,ra);elsez(a,b)=0;endendendxa=-1:h:1;ya=-1:h:1;[x,y]=meshgrid(xa,ya);subplot(2,2,1);surf(x,y,z);xlabel('x');ylabel('y');zlabel('u');title('数值解')p=sin(pi*(1-r^2)/2);fori=-1:h:1forj=-1:h:1ra=(i^2+j^2)^0.5;a=(i+1)/h+1;a=int8(a);b=(j+1)/h+1;b=int8(b);ifra1pp(a,b)=subs(p,r,ra);elsepp(a,b)=0;endendendsubplot(2,2,2);surf(x,y,pp);xlabel('x');ylabel('y');zlabel('p');title('精确解')fori=1:1:2/h+1forj=1:1:2/h+1e(i,j)=abs(z(i,j)-pp(i,j));endendsubplot(2,2,3);surf(x,y,e);xlabel('x');ylabel('y');zlabel('e');title('误差')error=max(max(e));end