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

例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:100
u(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)
end
meshz(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列的矩阵
figure
subplot(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; %初温度为零
figure
subplot(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:200
uu(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; %初温度为零
figure
subplot(1,2,1) %初始状态
plot(x,uu(:,1),'linewidth',1);
axis([0,1,0,0.25]);
subplot(1,2,2) %演化图
for j=2:200
uu(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:time
C(:,t+1)=4i*(A\B(:,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:time
C(:,t+1)=4i*(A\B(:,t));
B(:,t+1)=C(:,t+1)-B(:,t);
end
mo=moviein(time); %实时动画播放
for j=1:time
plot(x,abs(B(:,j)).^2/norm(B(:,j)).^2,x,v/2);
mo(:,j)=getframe;
end
movie(mo)

%例5：（书上例题1）两端固定的弦，初速为零，初位移不为零。
clear all
a=2;dt=0.05;dx=0.1;c=4*dt^2/dx^2;
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 all
N=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));
hplot(x,u(:,1),'linewidth',3);
axis([0,1,-0.05,0.05]);
set(h,'EraseMode','xor','MarkerSize',18)
for k=2:N
set(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 jxj
clear;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)];
end

function wtx=wfun(N,t)
x=0:0.001:1; a=1; wtx=0;
for I=1:N
if I~=7
wtx=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);
else
wtx=wtx+0.05/7*cos(I*pi*a*t).*sin(I*pi*x);
end
end

%例8：显示差分公式（迭代法）解平面温度场
u=zeros(100,100); u(100,)=10;
uold=u+2.5; unew=u;
for k=1:100
if 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;
end
end
surf(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:100
for i=2:19
for j=2:29
ph=(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);
end
end
end
colormap([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:10
Z2=5^4*5*((-1)^n*n^2*pi^2+2-2*(-1)^n)*...
sinh(n*pi.*Y/5).*sin(n*pi.*X/5)/...
(n^5*pi^5*sinh(n*pi*5/10));
Z1=Z1+Z2;
end
Z=Z1+X.*Y.*(5^3-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/2^2./T);
js=trapz(F,3);
waterfall(X(:,:,1), T(:,:,1), js)
figure
h=plot(xx',js(1,:))
set(h,'erasemode','xor');
for j=2:10
set(h,'ydata',js(j,:));
drawnow;
pause(0.1)
end

%例12：有限长细杆在第一类边界条件下的热传导问题
function jxj
N=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;
end
function u=rcdf(N,t)
x=0:0.21:20; u=0;
for k=1:2*N
cht=2/k/pi*(cos(k*pi*10/20)-...
cos(k*pi*11/20))*sin(k*pi*x./20);
u=u+cht*exp(-k^2*pi^2*10^2/400*t)); %a=10,l=20
end

%例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:60
lu(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 on
for j=2:length(t)
pause(0.01)
set(h,'ydata',jf(j,:));
drawnow;
end

%例15：半径为r0的球面径向速度分布已知，求球面所发射的稳恒声振动的速度势u
clear
r0=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.03
u=real(v0/2*r0^3*(-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 square
pause(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 3^0','P 3^1','P 3^2','P 3^3');

%例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)<0
h=interp1(y(k:k+1), x(k:k+1),0);
LD=[LD,h]
end
end

%例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