2u 1 u 1 2u 2 2 0 u 2 ku Hu f ; 0 q sin 0 其中： f 其他; 0
该问题的解析解为 ：
q 1 q 2q 2n u sin 2 n1 cos 2n 2 2nk Ha 1 4n H k Ha 2 n 1 a
先求其本征解。设函数 u(x,t) 具有变量分离形式：
u ( x, y) X ( x)Y ( y)，则上述方程可以写为：
X X 0; X 0 0; X a sin 3πy / b Y Y 0; Y 0 0; Y b / 2 * sin 2πx / a sin 4πx / a
考虑到对称性，解析解与φ 无关，因此在二维平面上画出等势线， 由于解为分段函数，因此分别算出并用不同颜色作图。
%ex404; (p100) 圆盘引力势； clear; a=0.35; GM=1/4; ri=0:1/200:a; ro=a:1/200:1; th=(0:0.01:2)*pi; z=cos(th); ui=0; uo=0; for k=0:2:20 fun=legendre(k,z); f=fun(1,:); f0=f(1,51); Ri=2/a^2*(1/(k-1)+1/(k+2)).*ri-2/a/(k-1)*(ri/a).^k; Ro=1/a/(k/2+1)*(a./ro).^(k+1); [R,PH]=meshgrid(Ri,f); ui =ui-GM*f0*R.*PH; [R,PH]=meshgrid(Ro,f); uo=uo-GM*f0*R.*PH; end; [ro,T]=meshgrid(ro,th);Yo=ro.*cos(T);Xo=ro.*sin(T); [ri,T]=meshgrid(ri,th); Yi=ri.*cos(T); Xi=ri.*sin(T); figure(1);contour(Xo,Yo,uo,'b');hold on; axis equal; contour(Xi,Yi,ui,'r'); title('圆盘引力势等势线');
由于轴对称性，可取x=0，可用trapz求先出电势，再用 gradient(梯度)求出电场，最后用streamline画出电力线
%ex403(p97) clear; a=1; b=0.11; y=-4:b:4; z=y; phi=pi*(0:1/100:2); [Y,Z,PHI]=meshgrid(y,z,phi); r=sqrt((0-a*cos(PHI)).^2+(Y-a*sin(PHI)).^2+Z.^2); dV=1./r; V=trapz(dV,3); [Ey,Ez]=-gradient(V,0.5); figure(2);subplot(2,2,1); contour(Y(:,:,1),Z(:,:,1),V,10); subplot(2,2,3); [Sy,Sz]=meshgrid(-4:.2:4,[-0.1,.1]); box on; streamline(Y(:,:,1),Z(:,:,1),Ey,Ez,Sy,Sz); x=0:b:3; th=pi*(0:1/20:2); [X,Y,Z,TH]=ndgrid(x,x,x,th); r=sqrt((X-a*cos(TH)).^2+(Y-a*sin(TH)).^2+Z.^2); dV=1./r; V=trapz(dV,4); [Ex,Ey,Ez]=gradient(-V,0.5); [X,Y,Z]=meshgrid(x); [Sx,Sy,Sz]=meshgrid(0:.5:3,0:.5:3,0.1); subplot(2,2,2); streamslice(X,Y,Z,Ex,Ey,Ez,[],[],0.1); box on; axis([0,3,0,3]);xlabel('Z_0=0.1'); x=cos(th);y=sin(th);z=zeros(1,length(th)); subplot(2,2,4); plot3(x,y,z,'linewidth',3,'color','r');hold on; axis([-3,3,-3,3,-3,3]);h1=streamline(X,Y,Z,Ex,Ey,Ez,Sx,Sy,Sz); h2=copyobj(h1,gca); rotate(h2,[1,0,0],180,[0 0 0]);
其分布如下页左图所示，右图则为偏微分方程工具箱 算出的结果。
用偏微分方程工具箱求解Laplace方程的步骤
1、在matlab命令窗中键入：pdetool 2、在弹出界面中用第二行的左边5个按钮之一选定求解边界. 3、用第二行的左边第6个按钮设定边界条件. 4、用第二行的左边第7个按钮设定求解方程
2 sin 2 x / a sinh 2 y / a / sinh 2 b / a if 2 / a 2 2 u x, y XY sin 4x / a sinh 4y / a / sinh 4b / a if 4 / a 2 sin 3y / b sinh 3x / b / sinh 3a / b if 3 / b 2
第五章 拉普拉斯方程与泊松方程
5.1 二维拉普拉斯方程
u u 2 0 2 x y
2 2
5.2 三维拉普拉斯方程 5.3 泊松方程与格林函数
u 0
u f r
5.1
二维拉普拉斯方程
5.1.1 矩形区域的拉普拉斯方程
u xx u yy 0 (0 x a;0 y b) u x 0, y 0; u x a, y sin 3y / b ; u x, y 0 0; u x, y b sin 3x / a cosx / a ;
(如选椭圆形并设定系数）. Nhomakorabea5、用第二行的左边第8、9个按钮剖分求解区域网格. 6、用第二行的左边第11个按钮画图.
5.1.2 圆形区域拉普拉斯方程，阳光照射的圆柱
半径为a,表面熏黑的长圆柱体，在温度为零度的 空气中受到垂直于柱轴的阳光照射，热流的强度为q, 求柱内温度分布。 因温度分布式稳定的，该问题的定界问题为：
2. 解析解的可视化 该定解问题的解析解----电势分布为： 2l q !r l 2l 1 2l P2l cos 2 l!l! a a l 1 u r , 2 l 1 ! a q l 2l 1 2l P2l cos a 2 l!l! r l 1
r a r a
由于对称性， 解与角φ无关，该解析解可以分为圆 内、圆外两部分，两部分分别计算，再将两部分相 加。程序如下页；相应的等位线在后页。
%ex402(p94) clear; a=0.5; q=1; x=[-3*a:0.02:3*a]; [X,Y]=meshgrid(x); [theta,r]=cart2pol(Y,X); rout=r; rout(find(rout<a))=NaN; rin=r; rin(find(rin>a))=NaN; Uin=q/a; Uout=q./rout; rin=rin/a; rout=a./rout; for k=1:20 fun=legendre(2*k,cos(theta)); rfun=q/a*squeeze(fun(1,:,:)); ck=(-1)^k*prod(1:2*k)/2^(2*k)/(prod(1:k))^2; ukin=ck*rin.^(2*k); Uin=Uin+ukin.*rfun; ukout=ck*rout.^(2*k+1); Uout=Uout+ukout.*rfun; end; figure(1);contour(X,Y,Uout,20,'r');hold on; contour(X,Y,Uin,15,'b');title('等势线'); figure(2);surf(X,Y,Uout);hold on;surf(X,Y,Uin);
u x, y

2
sin 2x / a sinh 2y / a / sinh 2b / a sin 4x / a sinh 4y / a / sinh 4b / a

2
sin 3y / b sinh 3x / b / sinh 3a / b
3. 用偏微分方程工具箱求解
5.2.3 均匀圆盘的引力势
均匀圆盘的半径为a，质量为M，求它在周围空间中的引力势。 定解问题：u 4GQr , 泊松方程
u 0、u 有解 u r 0有解、u r 0 其中： Qr , M / 2 H r a ; a 2 r 1 0 r a H r a 0 r a


求上述温度分布的程序如下，相应的分度分布如下页 上面两图，用偏微分方程工具箱所得分布可作为对比
%ex401(p91) clear;a=1; H=1.5; k=0.3; h=0.2; q=5; N=20;r=0:0.05:1; phi=0:pi/30:2*pi; [TH,R]=meshgrid(phi,r); %构造网格 [X,Y]=pol2cart(TH,R); u=q/H/pi+1/(k+H*a)*q/2*R.*sin(TH); for n=1:N dd=2*q/pi/a^(2*n-1)/(2*n*k+a*H)/(1-4*n^2); zz=dd*R.^(2*n).*cos(2*n*TH); u=u+zz; end; figure(1);surfc(X,Y,u); figure(2);contour3(X,Y,u,20);