スペクトル法で時間発展拡散方程式(1D)

勉強のためやってみたので、Octaveのスクリプトを。
もちろん、間違ってても責任は持ちません。

function laplacian = calc_lap(nDim, dif_coef, cur_func)
   for i=1:1:nDim
      laplacian(i) = -1.0*dif_coef*i^2*cur_func(i);
   end
endfunction
nGrid=128;
nSteps=2000;
deltaT=0.0001;
dif_coef = 1.0;
#Define initial condition
initial_value=zeros(1,nGrid);
initial_value(1,nGrid/2)=1.0;
init_fft=fft(initial_value);
#Begin time integral loop
cur_coef=zeros(1,nGrid);
hlf_coef=zeros(1,nGrid);
nxt_coef=zeros(1,nGrid);
lap_coef=zeros(1,nGrid);
cur_coef = init_fft;
#hold on;
#plot(abs(initial_value));
for iStep=1:1:nSteps
# 2-stage Runge-Kutta time integral
   lap_coef = calc_lap(nGrid, dif_coef, cur_coef);
   for i = 1:1:nGrid
      hlf_coef(i) = cur_coef(i) + deltaT/2.0 * lap_coef(i);
   end
   lap_coef = calc_lap(nGrid, dif_coef, hlf_coef);
   for i = 1:1:nGrid
      nxt_coef(i) = cur_coef(i) + deltaT * lap_coef(i);
   end
   cur_coef= nxt_coef;
   if(mod(iStep,50) == 0)
     cur_val = ifft(cur_coef);
     hold on;
     plot(real(cur_val));
   endif
end
plot(real(cur_val));
hold off;
print -deps time_diffusion