波动方程初值问题MATLAB程序

bodongchuzhi.m
a=1;
xmin=0;
xmax=2*pi;
supt=2;
nx=32;
nt=100;
dx=(xmax-xmin)/nx;
dt=supt/nt;
x=xmin:dx:(xmax-dx);
x=x';
t=0:dt:supt;
u=zeros(nx,nt+1);
for n=1:nx
    u(n,1)=sin(x(n));
end
for n=2:nx-1
    u(n,2)=1/2*a^2*(dt/dx)^2*(sin(x(n-1))+sin(x(n+1)))+(1-a^2*(dt/dx)^2)*sin(x(n))+dt*cos(x(n));
end
u(1,2)=1/2*a^2*(dt/dx)^2*(sin(x(nx))+sin(x(2)))+(1-a^2*(dt/dx)^2)*sin(x(1))+dt*cos(x(1));
u(nx,2)=1/2*a^2*(dt/dx)^2*(sin(x(nx-1))+sin(x(1)))+(1-a^2*(dt/dx)^2)*sin(x(nx))+dt*cos(x(nx));
r=a^2*dt^2/dx^2;
b=r;
c=2*(1-r);
Dxx=spdiags(kron([b,c,b],ones(nx,1)),-1:1,nx,nx);
Dxx(1,nx)=b;
Dxx(nx,1)=b;
for n=2:nt
    u(:,n+1)=Dxx*u(:,n)-u(:,n-1);
end