Contents
% This code simulates a diffusion equation on a 1D spatial domain [-L,L]; % It is a {\bf Fickian} diffusion % A reaction term F(u) can be added. For instance F(u)=ru(1-u/K). % Determine a function u(t,x) as a matrix where lines correspond to % different time steps and columns correspond to different space locations. % Initial Conditions are provided in the variable CI and plotted on % Figure(1) % Boundary conditions are Dirichlet conditions at infinity (u->0 when % x->+/-\infty % The diffusion coefficient is space dependent and given in function % DiffCoef.m, its derivative is needed and given in dDiffCoef.m % Jean-Christophe POGGIALE - Mars 2016 clear all;clc; global A B M L p reac %Parameter setting A=0.2;B=0.2; % Define the diffusion coefficient tmax=20; % Time length of the simulation dt=0.1; % Time step dx=0.05; % Space step L=10; % Domain=[-L,L] x=-L:dx:L; % Define the space mesh r=2;K=2;p=[r;K]; % Parameters of the reaction term, inserted in a vector p % OPTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BC0flux=0; % if =1 then no flux at the boundaries -L and L, else L and -L are not boundaries reac=0; % if =1 then the reaction term is included in the simulation SolAnalyt=1; % if =1, compares the numerical solution to the analytical solution %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Initial conditions x0=0;% Center of the initial distribution sigma=0.1; CI=exp(-((x-x0)/sigma).^2); % Gaussian curve around x_0 CI=CI/(sum(CI)*dx); % Normalization, the total amount of matter is 1 % Calculus of the diffusion matrix D=DiffCoef(x); dD=dDiffCoef(x); M=diag(-2*D/dx^2); MU=diag((dD(1:end-1)/2+D(1:end-1)/dx)/dx,1); ML=diag((D(2:end)/dx-dD(2:end)/2)/dx,-1); M=M+ML+MU; % Neuman Boundary conditions for a closed domain if BC0flux==1 M(1,2)=2*D(1,2)/dx^2; M(end,end-1)=2*D(end,end-1)/dx^2; end % Running simulation [t u]=ode45(@ModelDiff1D,0:dt:tmax,CI); % Analytical Solution t0=sigma^2/(2*DiffCoef(x0)); i=1;UA=[]; for i=0:length(t) uA=GaussFunct(t0+i*dt,x,x0,DiffCoef(x0),t0); UA=[UA;uA]; end
GRAPHS
figure(1) % Initial Condition plot(x,CI) xlabel('Space');ylabel('initial u distribution'); figure(2) % Diffusion coefficient as a function of space plot(x,D); ylim([0 max(D)]) xlabel('Space');ylabel('Diffusion coefficient'); figure(3) % dynamics of u on the spatial domaine for i=1:length(t) if SolAnalyt==1 plot(x,u(i,:),x,UA(i,:)); else plot(x,u(i,:)); end xlim([-L L]); ylim([0 2]) xlabel('Space');ylabel('u(t,x)') F(i)=getframe; end
![](Diffusion_01.png)
![](Diffusion_02.png)
![](Diffusion_03.png)
figure(4) % Total amount of matter versus time S=sum(u')*dx;S=S'; plot(t,S) ylim([0 max(S)+1]) xlabel('Time');ylabel('U_T'); title('Total amount of matter versus time')
![](Diffusion_04.png)