波动方程FDM，MATLAB

Question

我想实现MATLAB代码解决了波动方程，我的功能看起来是这样的：

function [x,t,w] = wave_eqn(xl,xr,yb,yt,M,N,f,l,r,p) 
% input: space interval [xl,xr], time interval [yb,yt] 
% number of space steps M, number of time steps N 
% output: solution w 
D=2;      % diffusion coefficient 
h=(xr-xl)/M; k=(yt-yb)/N; m=M-1; n=N; 
sigma=D*k/(h*h); 
a=diag(1-2*sigma*ones(m,1))+diag(sigma*ones(m-1,1),1); 
a=a+diag(sigma*ones(m-1,1),-1); % define matrix a 
lside=l(yb+(0:n)*k); rside =r(yb+(0:n)*k); 
w(:,1)=f(xl+(1:m)*h)'; % initial conditions 
for j=1:n 
    w(:,j+1)=a*w(:,j)-w(:,j-1)+sigma^2*[lside(j);zeros(m-2,1);rside(j)]; 
end 
w=[lside;w;rside];  % attach boundary conds 
x=(0:m+1)*h;t=(0:n)*k; 
% view(60,30);axis([xl xr yb yt -1 1]) 
end 

%source: numerical analysis 2nd edition 

我一直在该方程与W环形收到错误（：，j-1）表示：下标索引必须是真正的正整数或逻辑。

我不太清楚如何解决这个问题。还应该注意的是，f，p，l，r都是x和t的输入函数。我使用了一个用于编写热代码的模板，但我不确定如何实现第四个函数，p。谢谢。

Answer 1

您正在循环使用1:n，因此j-1=0和0不是Matlab中的有效索引。

取而代之的是，从2:n-1循环到j+1的术语。

你已经宣布你的初始条件w(:,1)，但你的数值方法需要2个前面的步骤，所以你还需要一个初始条件分配给w(:,2)，可能就这么简单w(:,2) = w(:,1)。然后循环使用：

for j=2:n-1 
    w(:,j+1)=a*w(:,j)-w(:,j-1)+sigma^2*[lside(j);zeros(m-2,1);rside(j)]; 
end