拉普拉斯方程与泊松方程课件.ppt

1、5.1 二维拉普拉斯方程二维拉普拉斯方程5.2 三维拉普拉斯方程三维拉普拉斯方程5.3 泊松方程与格林函数泊松方程与格林函数02222yuxu0u rfu先求其本征解。设函数 u(x,t)具有变量分离形式：，则上述方程可以写为：5.1 5.1 二维拉普拉斯方程二维拉普拉斯方程 5.1.1 矩形区域的拉普拉斯方程矩形区域的拉普拉斯方程;/cos/3sin,;00,;/3sin,;0,0)0;0(0axaxbyxuyxubyyaxuyxubyaxuuyyxx)()(),(yYxXyxu axaxbYYYYbyaXXXX/4sin/2sin*2/;00;0/3sin;00;0222/3 if /3s

2、inh/3sinh/3sin/4 if /4sinh/4sinh/4sin2/2 if /2sinh/2sinh/2sin2,bbabxbyaabayaxaabayaxXYyxu /3sinh/3sinh/3sin /4sinh/4sinh/4sin2 /2sinh/2sinh/2sin2,babxbyabayaxabayaxyxu其分布如下页左图所示，右图则为偏微分方程工具箱算出的结果。用偏微分方程工具箱求解用偏微分方程工具箱求解Laplace方程的步骤方程的步骤1、在matlab命令窗中键入：pdetool2、在弹出界面中用第二行的左边5个按钮之一选定求解边界选定求解边界.3、用第二行的左

3、边第6个按钮设定边界条件设定边界条件.4、用第二行的左边第7个按钮设定求解方程设定求解方程 (如选椭圆形并设定系数）如选椭圆形并设定系数）.5、用第二行的左边第8、9个按钮剖分求解区域网格剖分求解区域网格.6、用第二行的左边第11个按钮画图画图.5.1.2 圆形区域拉普拉斯方程，阳光照射的圆柱圆形区域拉普拉斯方程，阳光照射的圆柱;011022222fHukuuuuu其中：其中：;00 sin其他qf 半径为a,表面熏黑的长圆柱体，在温度为零度的空气中受到垂直于柱轴的阳光照射，热流的强度为q,求柱内温度分布。因温度分布式稳定的，该问题的定界问题为：该问题的解析解为该问题的解析解为 ：122122

4、cos4122sin21nnnnnHankaqqHakHqu%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*n2);zz=dd*R.(2*n).*cos(2*n*TH);u=u+zz;end;figure(1);surfc(X,Y,u);

5、figure(2);contour3(X,Y,u,20);求上述温度分布的程序如下，相应的分度分布如下页求上述温度分布的程序如下，相应的分度分布如下页上面两图，用偏微分方程工具箱所得分布可作为对比上面两图，用偏微分方程工具箱所得分布可作为对比5.1.3 云与大地之间的电缆云与大地之间的电缆带电的云与大地之间存在一个均匀的电场，平行与大带电的云与大地之间存在一个均匀的电场，平行与大地的电缆相当于一根无穷长导体。在平行于电场的方地的电缆相当于一根无穷长导体。在平行于电场的方向作垂直于电缆的截面，研究该截面上的电势分布。向作垂直于电缆的截面，研究该截面上的电势分布。该问题可以用如下方程加以描述：该问

6、题可以用如下方程加以描述：00circlesquarey;uuu选择偏微分方程工具箱求解，注意：求解区域是两区选择偏微分方程工具箱求解，注意：求解区域是两区域之差；区域内选择椭圆型方程。域之差；区域内选择椭圆型方程。5.2 5.2 三维拉普拉斯方程三维拉普拉斯方程 5.2.1 5.2.1 静电场中的介质球的电场分布静电场中的介质球的电场分布 在场强为在场强为E的均匀静电场中放置半径为的均匀静电场中放置半径为R0的均匀介的均匀介质球，球的介电常数质球，球的介电常数(电容率电容率)为为，求介质球内外的电，求介质球内外的电场强度分布。场强度分布。定界问题可以用下列方程加以描述：定界问题可以用下列方程

7、加以描述：2.02 ;0.y;u usquare；圆外圆内圆外1=2，圆内2=2 5.2.2 带有电荷的细圆环的电势分布带有电荷的细圆环的电势分布半径为半径为a 的均匀带电细圆环，电的均匀带电细圆环，电荷线密度为荷线密度为440/a,取无穷远处取无穷远处电势为电势为0,求空间任一点的电势求空间任一点的电势.1.直接积分法直接积分法20222200sincos 41,zayaxdradzyxU由于轴对称性，可取x=0，可用trapz求先出电势，再用gradient(梯度)求出电场，最后用streamline画出电力线%ex403(p97)clear;a=1;b=0.11;y=-4:b:4;z=y

8、;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

9、*(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,leng

10、th(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);h3=copyobj(allchild(gca),gca);rotate(h3,0,1,0,180,0 0 0);2.解析解的可视化解析解的可视化该定解问题的解析解该定解问题的解析解-电势分布为：电势分布为：arPralllaqarParlllaqrulllll

11、lllll cos!2!21 cos!2!21,212122212由于对称性，由于对称性，解与角解与角无关，无关，该解析解可以分为圆该解析解可以分为圆内、圆外两部分，两部分分别计算，再将两部分相内、圆外两部分，两部分分别计算，再将两部分相加。程序如下页；相应的等位线在后页。加。程序如下页；相应的等位线在后页。%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(routa)=NaN;Uin=q/a;Uout=q./rout;rin=rin/a;rout=

12、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

13、);hold on;surf(X,Y,Uin);3.用偏微分方程工具箱求解用偏微分方程工具箱求解 5.2.3 均匀圆盘的引力势均匀圆盘的引力势均匀圆盘的半径为均匀圆盘的半径为a，质量为，质量为M，求它在周围空间中的引力势。，求它在周围空间中的引力势。定解问题：定解问题：其中：其中：下面分别用两种方法求解：下面分别用两种方法求解：1.解析法绘图，解析法绘图，2.PDE工具箱求解工具箱求解000,4ruruuurGQu有解、有解、00 1 ;2/,2arararHarHraMrQ泊松方程1.解析解的可视化解析解的可视化问题的解析解为：考虑到对称性，解析解与无关，因此在二维平面上画出等势线，由于解为

14、分段函数，因此分别算出并用不同颜色作图。arPParlarllaGMarPPralaGMrullllllll if cos01212211212 if cos011,022202212%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/a2*(1/(k-1)+1/(k+2).*ri-2/a/(k-1)*(ri/a).k

15、;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(圆盘引力势等势线圆盘引力势等势线)

16、;2.用用PDE工具箱求解工具箱求解 选一长为选一长为3，宽为，宽为2的矩形，内画一小矩形表示板，的矩形，内画一小矩形表示板，长长(直径直径)取为取为0.8,厚度取为厚度取为0.1，中心位于原点。，中心位于原点。外边界条件是外边界条件是Dirichlet型：型：h=1,r=0;内边界条件自内边界条件自动衔接。动衔接。方程选为椭圆型，在圆板内参数取为：方程选为椭圆型，在圆板内参数取为：c=1,a=0,f=1/0.352./sqrt(x.2+y.2);圆板外参数取为：圆板外参数取为：c=1,a=0,f=0。5.2.4 环形电流的磁感应强度环形电流的磁感应强度 半径为a，通有电流I的圆环，求其在空间

17、任一点所产生的磁感应强度B.为了画出磁力线，需用到画场线的专用命令(streamline)。1.解析解解析解 arPParlrIarPPralrIrBarPParrIarPPrarIrBllllllllllllllllr if ;cos01212 if ;cos02212,if ;cos02 if ;cos02,01121121200112112220012112120012112220%ex406;(p104)%环形电流的磁感应强度环形电流的磁感应强度clear;a=0.35;R=1;X,Y=meshgrid(0:0.1:R);t,r=cart2pol(Y,X);ri=r;ri(find(r

18、ia)=NaN;ro=r;ro(find(ro0),j=j+1;end;%粗估零点粗估零点;q=fzero(yy,j);D=D,q;j=j+1;%精估一阶贝塞尔函数零点；精估一阶贝塞尔函数零点；endb=1.5;u0=1;h=0.5;rho,z=meshgrid(0:b/30:b,0:h/20:h);A=u0*b2/(2*h)*z;for n=1:N;aa=4*u0*besselj(2,D(n)/(D(n)/b)2;bb=sinh(D(n)*z./b)/sinh(D(n)*h/b);cc=besselj(0,D(n).*rho/b)/(besselj(0,D(n)2;A=A-aa*bb.*cc

19、;end;figure(1);contour(0:b/30:b,0:h/20:h,A,25);xlabel(radical rho);ylabel(high z);colorbar;5.2.6 柱体内温度场分布之二柱体内温度场分布之二(J的应用的应用)半径为 0高度为L的圆柱体，其侧面和下底温度保持0度，上底温度分布为 。求圆柱体内稳定的温度分布。定解问题为：其解析解为：其中 为零阶贝塞尔函数的第n个零点，J0和J1为0、2阶贝塞尔函数取0=0.2，L=0.5,用下列程序画图，结果如后页所示。;,;00,;,0 ;0,;020Lzuzuzuzuu有限2 010000000012020/sinh

20、/sinh412,nnnnnnnxJxxJLxzxxzu 0nx%ex411;(p117)%温度场分布之二温度场分布之二；clear;D=;N=50;j=1;yy=inline(besselj(0,xx),xx);for k=1:N;while(besselj(0,j)*besselj(0,j+1)0),j=j+1;end;%零点粗估零点粗估；q=fzero(yy,j);D=D,q;j=j+1;%精估精估零阶贝塞尔函数零点零阶贝塞尔函数零点；endr0=0.2;L=0.8;A=0;r,z=meshgrid(-r0:0.01:r0,0:0.01:L);for n=1:N;aa=2*r02/D(n

21、)/besselj(1,D(n)*(1-4/(D(n)2);bb=sinh(D(n)*z/r0)./sinh(D(n)*L/r0);cc=besselj(0,D(n)*r/r0);A=A+aa*bb.*cc;end;HH=10-5:3*10-4:10-3,.0015,.002:.001:.04;figure(1);contour(r,z,A,HH);xlabel(radical rho);ylabel(high z);colorbar;5.2.7 柱体内温度场分布之三柱体内温度场分布之三(I0 的应用的应用)半径为 0高度为L的圆柱体，侧面有均匀分布的热流流入，强度为q0,上、下底温度保持u0

22、。求圆柱体内稳定的温度分布。定解问题为：其解析解为：其中I0为第一类零阶修正（或第一类零阶虚宗量）贝塞尔函数。取u0=3;0=0.6,L=2.5,q0=1用下列程序画图，结果如后页所示;,;0,;,0 ;,;00000uLzuuzuzuqzuu有限LzlLlILlIlLLquzul/12sin/12/12124,000022200%ex412;(p120)%温度场分布之三温度场分布之三；clear;u0=3;r0=0.6;L=2.5;N=21;q0k=1;r,z=meshgrid(-r0:0.01:r0,0:0.01:L);A=u0;for n=1:2:N;b=n*pi/L;a=4*q0k/L

23、/b2;I1=(besseli(0,b*(r0-0.01)-besseli(0,b*r0)/(b*0.01);c=besseli(0,b*r)/I1;A=A+a*c.*sin(b*z);end;figure(1);contourf(r,z,A,20);xlabel(radical/rho);ylabel(high/z);colorbar;5.2.8 柱体内温度场分布之四柱体内温度场分布之四(I0 的应用的应用)半径为 0高度为L的圆柱体，上、下底温度分别保持u1、u2，侧面温度为:求圆柱体内稳定的温度分布。定解问题为：其解析解为：取u1=-5;u2=30;0=0.6,L=2.5,用下列程序画图

24、,结果如后页所示;,;0,;,0 ;,;0210uLzuuzuzuzfzuu有限LzpLpILpIpuzLuuuzup/12sin/12/121216,0000332121 zLLuLzLzuzf12222%ex413;(p122)%温度场分布之四温度场分布之四；clear;u1=-5;u2=30;L=2.5;r0=0.6;N=11;x=-r0:0.01:r0;y=0:0.005:L;r,z=meshgrid(x,y);u=u1+(u2-u1)/L*z;for n=1:2:N;a=16*u2/(n*pi)3;b=n*pi/L;u=u-a*besseli(0,b*r)./besseli(0,b*

25、r0).*sin(b*z);end;figure(1);contour(r,z,u,40);xlabel(radical/rho);ylabel(high/z);colorbar;5.2.9 电子透镜电子透镜(J0、I0 的应用的应用)如图所示，电子光学透镜的某一部件由两个中空圆柱筒组成，其电势分别为+v0和-v0，在圆柱中间隙缝的边缘处电势可近似为u=v0sin(z/2),求圆柱体内的电势分布。圆柱体两端边界条件可近似为u|=+-l=+-v0,圆柱体的半径为0。取柱坐标系，由于电势为z的奇函数，可以取通过缝隙中间的中点为柱坐标系的原点，则定解问题为：lvzvRuvlzuzvzuuzz0 z0

26、 ;2/sin;02sin0;000000令：u=v+w,将它分解为两个定解问题：其解析解分别为：取u1=-5;u2=30;0=0.6,L=2.5,用下列程序画图,结果如后页所示 lzvzzvvlzvzvv0 0 2/sin;0 ;00;0000lzplpIlpIpllplpvzvpp/sin/4/cos112,00012222100;00;000vvlzwzww 000000010100/sinh/sinh12,nnnnnnxJlxzxxJxvzw%ex414;(p125)%电子透镜；电子透镜；clear;a=;N=50;j=1;v0=100;yy=inline(besselj(0,xx),

27、xx);for k=1:N;while(besselj(0,j)*besselj(0,j+1)0),j=j+1;end;%粗估零点；粗估零点；q=fzero(yy,j);a=a,q;j=j+1;%精估贝塞尔函数零点；精估贝塞尔函数零点；end;h=1.5;r0=0.5;del=0.05;A=0;B=0;r,z=meshgrid(0:r0/20:r0,-h:h/20:h);for n=1:N;bb=besselj(0,a(n)*r/r0)/(a(n)*besselj(1,a(n);cc=sinh(a(n)*z/r0)/sinh(a(n)*h/r0);A=A+bb.*cc;K=n*pi/h;dd=

28、besseli(0,K*r)/besseli(0,K*r0);SS=h2-(2*n*del)2;ee=1/(n*pi)*(-1)(n+1)+h2/(SS+eps)*cos(K*del)*sin(K*z);B=B+dd.*ee;end;C=2*v0*(A+B);m=min(min(C);M=max(max(C);figure(1);surfc(r,z,C);axis(0,r0,-h,h,1.1*m,1.1*M);5.2.10 柱体外的电势柱体外的电势(K0 的应用的应用)半径为 0高度为L的导体圆柱壳，上下底和侧面被绝缘物质所隔离，柱壳侧面的电势为z，上下底接地。求柱体外的电势分布。定解问题为：

29、其解析解为：取 0=0.5,L=3,用下列程序画图,结果如后页所示;0 ;00;000LzuzuzuuLznLnKLnKnLzunn/sin/21,0000%ex415;(p128)%柱体外电势；柱体外电势；clear;a=;N=50;j=1;v0=100;h=3;r0=0.5;del=0.05;r,z=meshgrid(r0:r0/10:6*r0,0:h/60:h);A=0;for n=1:N;K=n*pi/h;bb=besselk(0,K*r)/besselk(1,K*r0).*sin(K*z);A=A+(-1)n*2/K*bb;end;figure(1);contour(z,r,A,20

30、);axis(0,h,r0,5*r0);ylabel(radical/r);xlabel(high/z);colorbar;5.3 5.3 泊松方程与格林函数泊松方程与格林函数 5.3.1 矩域的泊松方程矩域的泊松方程在矩形区域上在矩形区域上 求解求解:且在边界上的值为零，定界问题为：且在边界上的值为零，定界问题为：2/2/,0bybax,2yxu;02/;02/;0 ;00,2byubyuaxuxuyxu 15542233/sin/sinh2/sinh122112nnnaxnaynabnnbanxaxyu其解为：其解为：%ex416;(p130)%泊松方程之解泊松方程之解clear;a=5;

31、b=5;N=20;x=0:0.1:a;y=-b/2:0.1:b/2;X,Y=meshgrid(x,y);Z=X.*Y.*(a3-X.3)/12;for n=1:N;K=n*pi/a;aa=(-1)n*(n*pi)2+2*(1-(-1)n);bb=b/a*sin(K*X).*sinh(K*Y)/(K5*sinh(K*b/2);Z=Z+aa*bb;end;figure(1);mesh(X,Y,Z);colormap(hot);view(119,7);上图为解析解取前上图为解析解取前20项所得的结果。项所得的结果。下图是采用下图是采用pdetool所得的结果。边界条所得的结果。边界条件取为齐次件取为

32、齐次Dirichlet源项取为源项取为f=x.2.*y 5.3.2 球域的泊松方程与格林函数球域的泊松方程与格林函数 在半径为在半径为a的导体球内的导体球内(或球外或球外)，距球心为，距球心为r0处处放置电荷量为放置电荷量为40q的点电荷的点电荷,求它形成的静电场。求它形成的静电场。定界问题为：定界问题为：;0,4000arGzzyyxxqG20220202020cos2/cos2,rrraraqrarrrrqrQ由镜像法不难求得其解为：由镜像法不难求得其解为：%ex417;(p133)%球域格林函数；球域格林函数；clear;a=1;r0=2;r1=a2/r0;q=1;q1=-q*a/r0;

33、x,y=meshgrid(-3:0.01:3,-1:.01:4);theta,r=cart2pol(y,x);r(r=a)=NaN;u1=q./(sqrt(r02-2*r0*r.*cos(theta)+r.2)+eps);u2=q1./sqrt(r12-2*r1*r.*cos(theta)+r.2);u=u1+u2;t=pi*(0:1/20:2);sx=0.1*cos(t);sy=r0+0.1*sin(t);figure(1);contour(x,y,u,-1:20,3000);hold on;box on;colorbar;Ex,Ey=gradient(-u);streamline(x,y,

34、Ex,Ey,sx,sy);axis equal;plot(a*exp(i*t),r);axis(-a,a,-a,a);title(球域内点电荷的电场球域内点电荷的电场);figure(2);contour(x,y,u1,-1:20,3000);hold on;box on;colorbar;Ex,Ey=gradient(-u1);streamline(x,y,Ex,Ey,sx,sy);axis equal;plot(a*exp(i*t),r);axis(-a,a,-a,a);title(球域内源电荷的电场球域内源电荷的电场);figure(3);contour(x,y,u2,11);hold

35、on;box on;colorbar;Ex,Ey=gradient(-u2);streamline(x,y,Ex,Ey,sx,sy);axis equal;plot(a*exp(i*t),r);axis(-a,a,-a,a);title(球域内像电荷的电场球域内像电荷的电场);5.3.3 圆域的泊松方程与格林函数圆域的泊松方程与格林函数 在半径为在半径为a的导体圆柱外，距轴心为的导体圆柱外，距轴心为r0处放置单处放置单位长度电量为位长度电量为=0q的线电荷的线电荷,求它形成的静电场。求它形成的静电场。定界问题为：定界问题为：;0,00aGyyxxqG其解为：其解为：其中：其中：r0为源点位置，

36、为源点位置，r为场点位置为场点位置202202202020cos2lncos2ln2,rrraraarrrrraqrQ%ex419;(p136)%圆域格林函数；圆域格林函数；clear;a=1;r0=2;r1=a2/r0;q=1;x,y=meshgrid(-5:0.02:5,-5:.02:5);theta,r=cart2pol(y,x);r(r=a)=NaN;r(r=r0)=NaN;u1=q/(2*pi)*log(a./sqrt(r02-2*r0*r.*cos(theta)+r.2);u2=-q/(2*pi)*log(r1./sqrt(r12-2*r1*r.*cos(theta)+r.2);u

37、=u1+u2;t1=pi*(0:1/10:2);t=pi*(0:1/10:2);sx=0.1*cos(t);sy=r0+0.1*sin(t);figure(1);contour(x,y,u,15);hold on;box on;colorbar;Ex,Ey=gradient(-u);streamline(x,y,Ex,Ey,sx,sy);axis equal;plot(a*exp(i*t1),r);axis(-a,a,-a,a);title(圆域内线电荷的电场圆域内线电荷的电场);figure(2);contour(x,y,u1,15);hold on;box on;colorbar;Ex,Ey=gradient(-u1);streamline(x,y,Ex,Ey,sx,sy);axis equal;plot(a*exp(i*t1),r);axis(-a,a,-a,a);title(圆域内源电荷的电场圆域内源电荷的电场);figure(3);contour(x,y,u2,7);hold on;box on;colorbar;Ex,Ey=gradient(-u2);streamline(x,y,Ex,Ey,sx,sy);axis equal;plot(a*exp(i*t1),r);axis(-a,a,-a,a);title(圆域内像电荷的电场圆域内像电荷的电场);