计算物理--解偏微分方程编程问题.doc

Four short words sum up what has lifted most successful individuals above the crowd: a little bit more.-author-date计算物理-解偏微分方程编程问题计算物理-解偏微分方程编程问题%例1：有限长细杆的热传导的定解问题x=0:20; t=0:0.01:1; a2=10;r=a2*0.01;u=zeros(21,101);u(10:11,1)=1;for j=1:100u(2:20,j+1)=(1-2*r)*u(2:20,j)+r*( u(1:19,j)+ u(3:21,j);plot(u(:,j); axis(0 21 0 1); pause(0.1)endmeshz(u)%例2：非齐次方程的定解问题的解析解a2=50;b=5;L=1;x,t=meshgrid(0:0.01:1,0:0.000001:0.0005); Anfun=inline('2/L*(x-L/2).2.*exp(-b*x./2/a2).*sin(n*pi*x/L)','x','n','L','b','a2'); %定义内联函数u=0;for n=1:30 An=quad(Anfun,0,1,n,L,b,a2); %inline函数中定义x为向量，其它为标量 un=An*exp(-(n*n*pi*pi*a2/L/L+b*b/4/a2/a2).*t).*exp(b/2/a2.*x).*sin(n*pi*x/L); u=u+un; size(u);end mesh(x,t,u) %x,t，u都为501行101列的矩阵figuresubplot(2,1,1)plot(u(1,:)subplot(2,1,2)plot(u(end,:)%例3：非齐次方程的定解问题的数值解dx=0.01;dt=0.000001;a2=50;b=5;c=a2*dt/dx/dx;x=linspace(0,1,100)' %将变量设成列向量uu(1:100,1)=(x-0.5).2; %初温度为零figuresubplot(1,2,1) %初始状态plot(x,uu(:,1),'linewidth',1);axis(0,1,0,0.25);subplot(1,2,2) %演化图h=plot(x,uu(:,1),'linewidth',1);set(h,'EraseMode','xor')for j=2:200uu(2:99,2)=(1-2*c)*uu(2:99,1)+c*(uu(1:98,1)+ uu(3:100,1)-. b*dt/dx*(uu(3:100,1)-uu(2:99,1);uu(1,2)=0;uu(100,2)=0; %边界条件uu(:,1)=uu(:,2);uu(:,1)set(h,'YData',uu(:,1);drawnow;pause(0.01)end%此列还可以不用图形句柄，直接用plot画图dx=0.01;dt=0.000001;a2=50;b=5;c=a2*dt/dx/dx;x=linspace(0,1,100)' %将变量设成列向量uu(1:100,1)=(x-0.5).2; %初温度为零figuresubplot(1,2,1) %初始状态plot(x,uu(:,1),'linewidth',1);axis(0,1,0,0.25);subplot(1,2,2) %演化图for j=2:200uu(2:99,2)=(1-2*c)*uu(2:99,1)+c*(uu(1:98,1)+ uu(3:100,1)-. b*dt/dx*(uu(3:100,1)-uu(2:99,1);uu(1,2)=0;uu(100,2)=0; %边界条件uu(:,1)=uu(:,2);uu(:,1)plot(x,uu(:,1),'linewidth',1);axis(0,1,0,0.25);pause(0.01)end%例4：运用Crank-Nicolson公式解一维运动粒子贯穿势垒的薛定谔方程time=130;x=1:220; v(x)=0;v(105:115)=0.18; %势垒函数k0=0.5;d=10;x0=40; %波包的参数：平均动量，宽度和中心 B0=exp(k0*i*x).*exp(-(x-x0).2./(2*d.2); %初始的高斯波包B(:,1)=B0.'A=diag(-2+2i+v)+diag(ones(219,1),1)+diag(ones(219,1),-1);for t=1:timeC(:,t+1)=4i*(AB(:,t);B(:,t+1)=C(:,t+1)-B(:,t);plot(x,abs(B(:,t).2/norm(B(:,t).2,x,v/2); %画归一化后的概率密度和位势axis(0,300,0,0.1);pause(0.1)end%此列可采用实时动画播放time=130;x=1:220; v(x)=0;v(105:115)=0.18; %势垒函数k0=0.5;d=10;x0=40; %波包的参数：平均动量，宽度和中心 B0=exp(k0*i*x).*exp(-(x-x0).2./(2*d.2); %初始的高斯波包B(:,1)=B0.'A=diag(-2+2i+v)+diag(ones(219,1),1)+diag(ones(219,1),-1);for t=1:timeC(:,t+1)=4i*(AB(:,t);B(:,t+1)=C(:,t+1)-B(:,t);endmo=moviein(time); %实时动画播放for j=1:timeplot(x,abs(B(:,j).2/norm(B(:,j).2,x,v/2);mo(:,j)=getframe;endmovie(mo)%例5：（书上例题1）两端固定的弦，初速为零，初位移不为零。clear alla=2;dt=0.05;dx=0.1;c=4*dt2/dx2;x=0:dx:1; %弦长u1(1:11)=0; u2(1:11)=0; u3(1:11)=0; %预设三个矩阵 u1(2:6)=-x(2:6);u1(7:10)=-1.5+1.5*x(7:10); %初始条件plot(x,u1);axis(0,1,-1,1);pause(0.3)u2(2:10)=u1(2:10)+c/2*(u1(3:11)-2*u1(2:10)+u1(1:9); %初始速度plot(x,u2);axis(0,1,-1,1);pause(0.3)for k=2:20 %求解方程u3(2:10)=2*u2(2:10)-u1(2:10)+c*(u2(3:11)-2*u2(2:10)+u2(1:9);u1=u2; u2=u3;plot(x,u3);axis(0,1,-1,1);pause(0.3)end%例6：（书上例题2）两端固定的弦，初速为零，初位移不为零。clear allN=4010;dx=0.0024;dt=0.0005;c=dt*dt/dx/dx;x=linspace(0,1,420)' %弦长u(1:420,1)=0; %边界条件 u(181:240,1)=0.05*sin(pi*x(181:240)*7); %初始位移u(2:419,2)=u(2:419,1)+c/2*(u(3:420,1)-2*u(2:419,1)+u(1:418,1);h=plot(x,u(:,1),'linewidth',3);axis(0,1,-0.05,0.05);set(h,'EraseMode','xor','MarkerSize',18)for k=2:Nset(h,'XData',x,'YData',u(:,2) ;drawnow;u(2:419,3)=2*u(2:419,2)-u(2:419,1)+c*(u(3:420,2)-2*u(2:419,2)+u(1:418,2);u(2:419,1)=u(2:419,2); u(2:419,2)=u(2:419,3);end%例7：（书上例题3）弦自由振动的解析解function jxjclear;clc;N=50;t=0:0.005:2.0;x=0:0.001:1;ww=wfun(N,0);ymax=max(abs(ww);h= plot(x,ww,'linewidth',3);axis( 0, 1, -ymax, ymax)sy= ;for n=2:length(t)ww=wfun(N,t(n);set(h,'ydata',ww);drawnow;sy=sy,sum(ww);endfunction wtx=wfun(N,t)x=0:0.001:1; a=1; wtx=0;for I=1:Nif I=7wtx=wtx+0.05*( (sin(pi*(7-I)*4/7)-sin(pi*(7-I)*3/7). /(7-I)/pi-(sin(pi*(7+I)*4/7)-sin(pi*. (7+I)*3/7)/(7+I)/pi )*cos(I*pi*a*t).*sin(I*pi*x);elsewtx=wtx+0.05/7*cos(I*pi*a*t).*sin(I*pi*x);endend%例8：显示差分公式（迭代法）解平面温度场u=zeros(100,100); u(100,:)=10;uold=u+2.5; unew=u;for k=1:100if max(max(abs(u-uold)>=0.01 %取矩阵中最大的元素unew(2:99,2:99)=0.25*(u(3:100,2:99)+u(1:98,2:99)+u(2:99,3:100)+u(2:99,1:98);uold=u; u=unew;endendsurf(u)%例9：用松弛法解定解问题omega=1.5;x=linspace(0,3,30); y=linspace(0,2,20);phi(:,30)=sin(3*pi/2*y)'phi(20,:)=(sin(pi*x).*cos(pi/3*x);for N=1:100for i=2:19for j=2:29ph=(phi(i+1,j)+phi(i-1,j)+phi(i,j+1)+phi(i,j-1);phi(i,j)=(1-omega)*phi(i,j)+0.25*omega*(ph);endendendcolormap(0.5,0.5,0.5);surfc(phi)%例10：解析解与PDETOOL求解的比较。X,Y=meshgrid(0:0.1:5,-5/2:0.1:5/2);Z1=0;for n=1:1:10Z2=54*5*(-1)n*n2*pi2+2-2*(-1)n)*. sinh(n*pi.*Y/5).*sin(n*pi.*X/5)/. (n5*pi5*sinh(n*pi*5/10);Z1=Z1+Z2;endZ=Z1+X.*Y.*(53-X.3)/12;colormap(hot);mesh(X,Y,Z)view(119,7)%例11：无限长细杆的热传导问题xx=-10:.5:10;tt=0.01:0.1:1;tau=0:0.01:1;a=2;X,T,TAU=meshgrid(xx, tt, tau);F=1/2/2./sqrt(pi*T).*exp(-(X-TAU).2/4/22./T);js=trapz(F,3);waterfall(X(:,:,1), T(:,:,1), js)figureh=plot(xx',js(1,:)set(h,'erasemode','xor');for j=2:10set(h,'ydata',js(j,:);drawnow;pause(0.1)end%例12：有限长细杆在第一类边界条件下的热传导问题function jxjN=50; t=1e-5:0.00001:0.005; x=0:0.21:20;w=rcdf(N,t(1);h= plot(x,w,'linewidth',5);axis( 0, 20, 0, 1.5)for n=2:length(t)w=rcdf(N,t(n);set(h,'ydata',w);drawnow;endfunction u=rcdf(N,t)x=0:0.21:20; u=0;for k=1:2*Ncht=2/k/pi*(cos(k*pi*10/20)-.cos(k*pi*11/20)*sin(k*pi*x./20);u=u+cht*exp(-(k2*pi2*102/400*t); %a=10,l=20end%例13：（波动方程）无限长的自由振动的弦满足的振动方程 %初位移不为零,初速为零u(1:140)=0; x=linspace(0,1,140);u(61:80)=0.05*sin(pi*x(61:80)*7);uu=u;h=plot(x,u,'linewidth',3);axis(0,1,-0.05,0.05);set(h,'EraseMode','xor')for at=2:60lu(1:140)=0; ru(1:140)=0;lx=61:80-at; rx=61:80+at;lu(lx)=0.5*uu(61:80); ru(rx)=0.5*uu(61:80);u=lu+ru;set(h,'XData',x,'YData',u) ;drawnow;pause(0.1)end%例14：（波动方程）无限长的自由振动的弦满足的振动方程 %初位移为零,初速不为零t=0:0.005:8; x=-10:0.1:10; a=1;X,T=meshgrid(x,t);xpat=X+a*T;xpat(find(xpat<=0)=0;xpat(find(xpat>=1)=1;xmat=X-a*T;xmat(find(xmat<=0)=0;xmat(find(xmat>=1)=1;jf=1/2/a*(xpat-xmat);%jf=1/2/a*xpat; %波形向左移动%jf=1/2/a*xmat; %波形向右移动h=plot(x,jf(1,:),'linewidth',3) ; %画动画set(h,'erasemode','xor');axis(-10 10 -1 1)hold onfor j=2:length(t)pause(0.01)set(h,'ydata',jf(j,:);drawnow;end%例15：半径为r0的球面径向速度分布已知，求球面所发射的稳恒声振动的速度势uclearr0=0.2; v0=2;k=60; a=2;theta=linspace(0,2*pi,50);rho=0.2:0.1:4;Th,Rh=meshgrid(theta,rho);X,Y=pol2cart(Th,Rh);rh=sqrt(X.2+Y.2);th=atan(Y./X);for t=0:0.001:0.03u=real(v0/2*r03*(-1./rh.2+i*k./rh).*.cos(th).*exp(k*(rh+2*t)*i);surf(X,Y,u);view(-32,28)%contour(X,Y,u); axis(-4 4 -4 4); axis squarepause(0.5)end%例16：勒让德多项式的图像x=0:0.01:1;y1=legendre(1,x);y2=legendre(2,x);y3=legendre(3,x);y4=legendre(4,x);y5=legendre(5,x);y6=legendre(6,x);plot(x,y1(1,:),x,y2(1,:),x,y3(1,:),x,y4(1,:),x,y5(1,:),x,y6(1,:)%例17：勒让德函数图像（以x为变量）x=0:0.01:1;y=legendre(3,x);plot(x, y(1,:),'-',x,y(2,:),'-.',x,y(3,:),':',x, y(4,:),'-')legend('P 30','P 31','P 32','P 33');%例18：勒让德函数图像（以theta为变量）rho=legendre(1,cos(0:0.1:2*pi);t=0:0.1:2*pi;rho1=legendre(2,cos(0:0.1:2*pi);rho2=legendre(3,cos(0:0.1:2*pi);subplot(3,4,1), polar(t,rho(1,:)subplot(3,4,2), polar(t,rho(2,:)subplot(3,4,5), polar(t,rho1(2,:)subplot(3,4,6), polar(t,rho1(2,:)subplot(3,4,7), polar(t,rho1(2,:)subplot(3,4,9), polar(t,rho2(2,:)subplot(3,4,10), polar(t,rho2(2,:)subplot(3,4,11), polar(t,rho2(2,:)subplot(3,4,12), polar(t,rho2(2,:)%例19：贝塞尔(Bessel)函数y=besselj(0:3,(0:.2:10)');figure(1)plot(0:.2:10)',y)legend('J0','J1', 'J2','J3')%例20：内插法与贝塞尔函数的零点x=0:0.05:50; y=besselj(0,x); LD=;for k=1:1000,if y(k)*y(k+1)<0h=interp1(y(k:k+1), x(k:k+1),0);LD=LD,hendend%例21：诺伊曼(Neumann)函数y=bessely(0:1,(0:.2:10)');plot(0:.2:10)',y)grid on%例22：虚宗量贝塞尔函数I=besseli(0:1,(0.1:0.1:3)');plot(0.1:0.1:3)',I)%例23：虚宗量汉克尔函数K=besselk(0:1,(0.1:0.1:3)');plot(0.1:0.1:3)',k)%例24：球贝塞尔函数x=eps:0.2:15;y1=sqrt(pi/2./x).*besselj(1/2,x);y2=sqrt(pi/2./x).*besselj(3/2,x);y3=sqrt(pi/2./x).*besselj(5/2,x);y4=sqrt(pi/2./x).*besselj(7/2,x);plot(x,y1,x,y2,x,y3,x,y4)%例25：球诺伊曼函数x=0.8:0.2:15;y1=sqrt(pi/2./x).*bessely(1/2,x);y2=sqrt(pi/2./x).*bessely(3/2,x);y3=sqrt(pi/2./x).*bessely(5/2,x);y4=sqrt(pi/2./x).*bessely(7/2,x);plot(x,y1,x,y2,x,y3,x,y4)axis(1 10 -0.5 0.4)grid on-