MATLAB 偏微分方程式　数値シミュレーション

MATLABを用いた、2次元偏微分方程式(反応拡散方程式)の濃度変化のシミュレーションについて質問します。

MATLABで下記のプログラムを書き、シミュレーションを行いました。

ある偏微分方程式モデル(2変数)のシミュレーションを行うためのプログラムです。

function iplus = ip(a,N)
if a >= N
iplus = a - (N-1);
else
iplus = a + 1;
end

function iminus = im(a,N)
if(a <= 1)
iminus = a + (N-1);
else
iminus = a - 1;
end

%main

clear all;

%*各種変数の設定(時間刻み幅、格子間隔など)。

N = input('格子点の数を入力：');
L = 1.; % 系の辺
h = L/(N-1); % 格子間隔
time = input('時間刻み幅を入力：');
D = 10^(-7); % 拡散係数
coeff = D * time / h^2;
alpha = 3.0*10^(-2); % the late of release
beta = 1; % the stlength of cell responce
tau = 64*60; % the characteristic degradation time of chemoattractant

%* 初期条件と境界条件を設定

cell = input('細胞数を入力：');
nna = randn(N, N);
nnasigma = 0;

for i = 1:N
for j = 1:N
nnasigma = nnasigma + nna(i, j);
end
end

nn = zeros(N, N);
for i = 1:N
for j = 1:N
nn(i, j) = (nna(i, j) / nnasigma) * cell;
end 
end


nn_new = zeros(N, N);
vvx = zeros(N, N);
vvy = zeros(N, N);
vvx_new = zeros(N, N);
vvy_new = zeros(N, N);
cc = zeros(N, N);
cc_new = zeros(N, N);

xplot = (0:N-1) * h - L/2;
yplot = (0:N-1) * h - L/2;
iplot = 1;
nstep = input('時間刻みの数を入力：');
nplots = 50;
plot_step = nstep/nplots;

clf;

figure(1);
mesh(xplot,yplot,nn);
figure(1);
contourf(nn,1000,'LineColor','none')
colorbar;

%差分法 

for istep = 1:nstep
for i = 1:N
for j = 1:N
cc_new(i,j) = cc(i,j) + ...
coeff*(cc(ip(i,N),j) + cc(im(i,N),j) + cc(i,ip(j,N)) + cc(i,im(j,N)) -4*cc(i,j)) + ...
time*(alpha * nn(i,j) - (1/tau) * cc(i,j));
end
end

for i = 1:N
for j = 1:N
vvx_new(i,j) = vvx(i,j) - (time/(2*h)) * ((vvx(ip(i,N),j) - vvx(im(i,N),j)) * vvx(i,j) + ...
(vvy(ip(i,N),j) - vvy(im(i,N),j)) * vvy(i,j)) + ...
(time/(2*h)) * beta * (cc(ip(i,N),j) - cc(im(i,N),j));

vvy_new(i,j) = vvy(i,j) - (time/(2*h)) * ((vvy(i,ip(j,N)) - vvy(i,im(j,N))) * vvx(i,j) + ...
(vvy(i,ip(j,N)) - vvy(i,im(j,N)))*vvy(i,j)) + ...
(time/(2*h)) * beta * (cc(i,ip(j,N)) - cc(i,im(j,N)));
end
end

for i = 1:N
for j = 1:N
nn_new(i,j) = nn(i,j) + ...
(time/(2*h)) * (nn(ip(i,N),j) * vvx(ip(i,N),j) - nn(im(i,N),j) * vvx(im(i,N),j) + ...
nn(i,ip(j,N)) * vvy(i,ip(j,N)) - nn(i,im(j,N)) * vvy(i,im(j,N)));
end
end
nnp = nn;
nn = nn_new;

vvx = vvx_new;

vvy = vvy_new;

cc = cc_new;

csvwrite('gambadensity',nn);

sigmadens = 0;
sigmacount = 0;

if(rem(istep,plot_step) < 1)
nnplot = nn;
iplot = iplot + 1;
if(rem(iplot, 10) == 0)
for i = 1:N
for j = 1:N
sigmadens = nn(i,j) + sigmadens;
end
end
sigmacount = sigmadens * L^2;
sigmacount
figure(iplot);
mesh(xplot,yplot,nnplot);
xlabel('x'); ylabel('y'); zlabel('N(x,y)');
grid on

figure(iplot*10);
contourf(nnplot,1000,'LineColor','none');
colorbar;
end 
end
end