Волновое уравнение с 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

Я продолжаю получать ошибку в уравнении в цикле for с термином w(:,j-1), который говорит: индексы индекса должны быть либо действительными положительными целыми числами, либо логическими.

Я не совсем уверен, как решить эту проблему. Также следует отметить, что f, p, l, r - все введенные функции x и t. Я использовал шаблон для уравнения теплопроводности для создания этого кода, но я не уверен, как реализовать четвертую функцию, с. Благодарю.

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