偏微分方程学习

1、椭圆形方程差分方法应用带微分余项的Taylor 公式可得()()()()()()()()()()()()6665544'''3''2'7201202462h O c u h c u h c u h c u h c u h c hu c u h c u +++++++=+()()()()()()()()()()()()6665544'''3''2'7201202462h O c u h c u h c u h c u h c u h c hu c u h c u ++-+-+-=-将上两式相加，得()()()[]()()()()()()66644''23601221h O c u h c u h c u h c u c u h c u h +++=-+-+ 使用标准的二阶中心差分格式，可得2,1,,1,2x u u u u ji j i j i ji x ∆+-=-+δ()6664436012x O u x u x u x x xx ∆+∆+∆+=2,1,,1,2yu u u u ji j i j i ji y ∆+-=-+δ()6664436012y O u y u y u y y yy ∆+∆+∆+=我们可以分别应用符号四阶紧致逼近算子到方程1中的二阶导xx u 和yy u 。

离散的二维泊松方程可以表示为如下形式()6424121222122121121∆+∆+∆+=⎪⎪⎭⎫ ⎝⎛∆++⎪⎪⎭⎫ ⎝⎛∆+--O y x f u y u x y y x x ττδδδδ 1τ和2τ用来表示在理查德森外推过程中忽略掉的误差，()6∆O 表示六阶精度()66y x O ∆+∆。

应用符号算子，设等于0，略去高阶项，方程可以写成()()()42222422222222221211121121121121∆+⎥⎦⎤⎢⎣⎡∆+∆+=∆+⎪⎪⎭⎫ ⎝⎛∆+⎪⎪⎭⎫ ⎝⎛∆+=⎪⎪⎭⎫ ⎝⎛∆++⎪⎪⎭⎫ ⎝⎛∆+O f y x O f y x u x u y y x y x y x x y δδδδδδδδ如果我们设置网格比为y x ∆∆=/λ，则方程可以写成()()()()1,1,,1,2,,1,11,1,1,11,11,11,182-++-+-+--+--++++++∆=-+++++++j i j i j i j i j i j i j i j i j i j i j i j i j i f f f f x au u u b u u c u u u u d系数是()2110λ+=a ， 25λ-=b ， 152-=λc ， ()2/12λ+=d .当网格比为1时，差分格式为如下形式()()()()1,1,,1,2,,1,11,1,1,11,11,11,1822044-++-+-+--+--++++++∆=-+++++++j i j i j i j i ji j i j i j i j i j i j i j i j i f f f f x u u u u u u u u u用矩阵可表示为⎪⎪⎪⎭⎫ ⎝⎛1- 4- 1-4- 20 4-1- 4- 1-由于该差分格式是九点紧致格式，四阶收敛，所以误差可以表示为)()(4h O h U A e h =-= )16()2(42h O h U A e h =-=这里h e 表示在步长为h 的情况下的误差，A 为精确解，)(h U 表示在步长为h 的情况下的数值解。

由此可以看出，当步长降低一倍时，误差比值为16，这个可以有效的检验计算结果的准确性。

数值算例()(),10,10 y sin 12≤≤≤≤-=∆-y x e u x ππ()(),sin ,0y y u π=()(),sin ,1y e y u π=,10≤≤y(),00,=x u (),01,=x u .10≤≤x该问题的精确解为()y e u xπsin =.表1给出了不同步长取值下的误差，误差比值及计算时间。

表1 计算误差表误差比值接近16，因此可以看出计算结果是比较准确的。

图1、图2、图3分别为步长0.05时的精确解曲面，数值解曲面和误差曲面。

图1 精确解曲面图2 步长为0.05时的数值解曲面图3 步长为0.05时的误差曲面MATLAB程序function [ error ] = jiudian( h )%UNTITLED Summary of this function goes here% Detailed explanation goes heretic;n = 1/h;for i=1:n+1 %boundary condition u(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%coefficient matrix AB = sparse(20*eye(n-1));D = sparse(-4*eye(n-1));for i = 1:n-2B(i,i+1) = -4;B(i+1,i) = -4;endfor i = 1:n-2D(i,i+1) = -1;D(i+1,i) = -1;enda = cell(n-1);for i = 1:n-1for j = 1:n-1a{i,j} = sparse(zeros(n-1));endendfor i = 1:n-1a(i,i) = {B};endfor i = 1:n-2a{i,i+1} = D;a{i+1,i} = D;end%Cell Array to MatrixA = cell2mat(a);A = sparse(A);% right vector Kfor i = 1:n+1for j = 1:n+1f(i,j) = (pi*pi-1)*exp((i-1)*h)*sin(pi*(j-1)*h);endendfor i = 2:nfor j = 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);for i = 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);%solve XX = A\K;for i = 2:nfor j = 2:nu(j,i) = X((n-1)*(i-2) + (j-1));endend%plot X, ex = 0:h:1;y = 0:h:1;subplot(2,2,1); surf(x,y,u);xlabel('x');ylabel('y');zlabel('u');title('差分解')for i = 1:n+1for j = 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('精确解');for i = 1:n+1for j = 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 方法求边值问题()(),21sin21cos2222222r r r u -+-=∆-ππππ()Ω∈y x ,，,0=u (),,Γ∈y x其中(){}1|,22<+=Ωy x y x ，222y x r +=。

该边值问题的精确解为()()21sin,2r y x u -=π。

由于Ω关于原点是对称的，方程右端函数仅依赖于2r 。

在()Ω10H 上取一簇奇函数()()(),,12-=j j r r w y x φn j ,...,2,1=.为使01==r jφ，取(),12r r w -=因而()()(),1,122--=j j r r y x φn j ,...,2,1=.以n φφφ,...,,21为基底张成n 维子空间n V ，n V 中任一元素可表示为()()()().1,,11221∑∑=-=-==nj j j nj j j n r c ry x c y x u φ计算得()(),1,'⎪⎪⎭⎫ ⎝⎛=∇y x r r y x i i ϕφ()()()()()(),,,''''r r y x y x j i j i φφφφ⋅=∇⋅∇ ()()()()(),2,,,1dr drd dr d rdxdy y x y x a ji j i j i φφπφφφφ⎰⎰⎰=∇⋅∇=Ω,,1n j i ≤≤()()()()()()()⎰⎰⎰⎥⎦⎤⎢⎣⎡-+--==-Ω12222122,21sin 21cos 212,,,dr r r r rr dxdyy x f y x f ri i i πππππφφn i ≤≤1， 将上述两式代入到Ritz-Galerkin 方程组，计算即可。

表2给出了取不同n 值下的误差。

图4为精确解曲面，图5和图6分别为n=4时的数值解和误差曲面，图7和图8分别为n=8时的数值解和误差曲面。

表2 不同n值的误差图4 n=4精确解曲面图5 n=4的数值解曲面图6 n=4的误差曲面图7 N=8的数值解曲面图8 N=8的误差曲面MATLAB程序function [ error ] = yxyty( n )syms r;for i = 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); for i = 1:nfor j = 1:na(i,j) = 2*pi*int(r*d(i)*d(j),0,1);endendfor i = 1:nb(i) = 2*pi*int(r*p(i)*f,0,1);endb = b';c = a\b;c = double(c);u1 = 0;for i = 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;for i = -1:h:1for j = -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);if ra<1z(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);for i = -1:h:1for j = -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);if ra<1pp(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('精确解')for i = 1:1:2/h+1for j = 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。