扩散方程
一阶微分
二阶微分
引入标记
得到离散的扩散公式
整理可得
编程：
for i = 1:numTimeSteps
for j = 1:numSpaceSteps
Anew(i,j)=A(i,j)+c*(A(i,j+1)+A(i,j-1)-2*A(i,j));
end
end
重要：，如果，解不稳定。显式稳定性差。
% diffuse1.m - Program to solve the diffusion equation
% using the FTCS method.
clear; % clear memory
% help diffuse1.m % print help information
%%%%%%%% INITIALISATION %%%%%%%%%%
% set up input variables
tmax = 1; % window size in t
nt = 100; % number of grid points in t
delta_t = tmax/nt; % step-size in t
nplots = 21; % number of plot points in t
nplotstep = nt/nplots; % plot step-size in t
xmax = 2; % window size in x
nx = 11; % number of grid points in x
delta_x = xmax/(nx + 1); % step-size in x
kappa = 1; % coupling constant
c = kappa*delta_t/delta_x^2; % constant of propagation
% initial conditions and boundary conditions
A = zeros(nx,1); % allocates and zeros a column vector
A(floor((nx+1)/2)) = 1/delta_x; % sets up an initial delta 'spike' in A
iplot = 1; % initialise the plotting index
% set up plotting matrices
xplot = (1:nx)*delta_x - xmax/2; % row vector of x-coordinates
tplot(iplot) = 0; % assign first time-coordinate
Aplot(:,iplot) = A(:); % assign first value of A
%%%%%%%% MAIN LOOP %%%%%%%%%%
% main loop
for i = 1:nt
% evolve our solution using the FTCS method
A(2:(end-1)) = A(2:(end-1)) + c*( A(1:(end-2)) + A(3:end) - 2*A(2:(end-1)) );
% save current data for plotting
if (rem(i,nplotstep) < 1)
iplot = iplot + 1; % increment plotting index
Aplot(:,iplot) = A(:); % record current value of A for plotting
tplot(iplot) = i*delta_t; % record current time value
% print some progress information
% fprintf('%g out of %g steps completed\n',i,nt);
end
end
%%%%%%%% PLOTTING/OUTPUT %%%%%%%%%%
% plot stuff
figure(1) % opens a new figure window
subplot(121) % divides window into a 1x2 matrix
% and assigns next plot to the
% left-hand frame (frame 1)
mesh(tplot,xplot,Aplot) % wire mesh plot of the solution
view([-70 30]) % change the viewing direction
matlab求导xlabel('t') % label the x-axis
ylabel('x') % label the y-axis
title('Diffusion of a delta spike')% give the thing a title
subplot(122) % assigns next plot to right-hand
% frame (frame 2)
% find contours of 3D plot and assign to the vector cs
cs = contour(xplot,tplot,flipud(rot90(Aplot)),0:0.1:2.5);
xlabel('x') % label the x-axis
ylabel('t') % label the y-axis
clabel(cs,0:5) % label the contours
title('Contour of a delta spike') % give the thing a title
subplot(111) % return plotting mode to one
% figure per window
离散的扩散方程可重写为
矩阵是微分算子的离散形式，，是Kronecker delta functions。
由于是列向量的单元，可将扩散方程写为
为单位矩阵
；
% diffuseMatrix.m - Program to solve the diffusion equation
% using the FTCS method with a matrix representation.
clear; % clear memory
% help diffuseMatrix.m % print help information
%%%%%%%% INITIALISATION %%%%%%%%%%
% set up input variables
tmax = 1; % window size in t
nt = 100; % number of grid points in t
delta_t = tmax/nt; % step-size in t
nplots = 21; % number of plot points in t
nplotstep = nt/nplots; % plot step-size in t
xmax = 2; % window size in x
nx = 11; % number of grid points in x
delta_x = xmax/(nx + 1); % step-size in x
kappa = 1; % coupling constant
c = kappa*delta_t/delta_x^2; % constant of propagation
发布评论