Schnakenberg system
clear all ;
close all;
D1=1;
D2=.1;
N=200;
L=100;
c=0;
t_end=100;
x =linspace (c,L,N);
dx= x(2) -x(1);
dt=.001
t =1:dt:t_end;
u = zeros(length(x),length(t));
v = zeros(length(x),length(t));
u(:,1) = 0;
u(70:140,1) = 150;
v(:,1) =0;
v(70:140,1) = 150;
for n=1:length(t)-1
for i=2:length(x)-1
u(i,n+1) = u(i,n) + D1*dt/(dx^2)*(u(i-1,n)-2*u(i,n)+u(i+1,n));%+dt*a-dt*T_mat1(i,n)+dt*(T_mat1(i,n)^2*(T_mat2(i,n)));
v(i,n+1) = v(i,n) + D2*dt/(dx^2)*(v(i-1,n)-2*v(i,n)+v(i+1,n));%+dt*b-dt*(T_mat1(i,n)^2*(T_mat2(i,n)));
end
u(1,n+1) = u(1,n) +dt/2*(dx)*(u(2,n)-u(1,n));
v(1,n+1) = v(1,n) +dt/2*(dx)*(v(2,n)-v(1,n));
u(end,n+1)= u(end,n)+dt/2*(dx)*(u(end-1,n)-u(end,n));
v(end,n+1)= v(end,n)+dt/2*(dx)*(v(end-1,n)-u(end,n));
end
for n=1 :length(t)
Total_u(n)= sum(u(1:length(x)-1,n))*dx;
end
for n=1 :length(t)
Total_v(n)= sum(v(1:length(x)-1,n))*dx;
end
figure
[tt,xx] = meshgrid (t,x);
zz=tt.*xx;
mesh(xx,tt,u);
xlabel('X coordinate (m)')
ylabel('Time (sec)')
zlabel('Temperature (F)')
zlim([0,200])
figure
[tt,xx] = meshgrid (t,x);
zz=tt.*xx;
mesh(xx,tt,v);
xlabel('X coordinate (m)')
ylabel('Time (sec)')
zlabel('Temperature (F)')
zlim([0,200])
figure
hold on
plot(x, u(:,1),x,v(:,1)) % initial condition T(x,0) vs x
% plot(x_vec, T_mat(:, 500)) % T(x, t_end/3) vs x
%plot(x_vec, T_mat(:, round(length(t_vec)/4)))
plot(x, u(:, round(length(t)/3)),x,v(:, round(length(t)/3))) % T(x, t_end/3) vs x
plot(x, u(:, round(length(t)*(2/3))),x, v(:, round(length(t)*(2/3)))) % T(x, t_end*2/3) vs x
%plot(x_vec, T_mat(:, round(length(t_vec)*(3/4))))
plot(x, u(:, length(t)),x,v(:, length(t))) % final state T(x, t_end) vs x
legend('t=0','t=(1/3)t_{end}', 't=(2/3)t_{end}','t=t_{end}')
hold off