波动方程的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)项的for循环中的等式中得到错误:下标索引必须是实数正整数或逻辑索引。

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

Answer 1

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

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

您已经将初始条件声明为w(:,1)，但是您的数值方法需要前面的两个步骤，因此您还需要为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