波動方程式、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 

私はワットとループの中の方程式でエラーを得続けます（：、j-1）という言葉：添え字インデックスは、正の整数または論理でなければなりません。

この問題を解決する方法がわかりません。また、f、p、l、rはすべてxとtの入力関数であることにも注意してください。このコードを作成するための熱方程式のテンプレートを使用しましたが、第4の関数pの実装方法がわかりません。ありがとう。

Answer 1

1:nをループしていますので、j-1=0と0はMatlabの有効なインデックスではありません。

2:n-1のループは、j+1という用語にもなります。

あなたはすでにあなたの初期条件w(:,1)は宣言しているが、あなたはまた、おそらくw(:,2) = w(:,1)ような単純な、w(:,2)に初期条件を割り当てる必要がありますので、あなたの数値的方法は、2つの前の手順が必要です。ループを使用してループします。

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