Hi DO,
It is my understanding that the question is about solving the 1D diffusion equation using the Crank-Nicolson method
To do the same, you can please try the following process:
- Setup the spatial domain and grid where the spatial grid ‘x’ is created with a spacing ‘dx’.
- Once the same is created, the time step size ‘dt’ is chosen based on the stability condition of the FTCS (Forward Time Centred Space scheme) and ‘alpha’ is calculated which represents the weight of the diffusion term.
- Construct a matrix (A) for the implicit Crank-Nicolson scheme and a gaussian distribution is set as the initial condition. The solution is then advanced in time using an implicit scheme and the results are plotted every 500 time steps.
% Define the x-domain and x-grid
diffusionCoefficient = 1;
domainLength = 1;
numIntervals = 500;
numGridPoints = numIntervals + 1;
dx = 2 * domainLength / numIntervals;
x = -domainLength + (0:numIntervals) * dx;
% Time step parameters
numTimeSteps = 10000;
plotFrequency = 500;
dt = dx^2 / (2 * diffusionCoefficient);
alpha = dt * diffusionCoefficient / dx^2; % Crank-Nicolson parameter
% Construct the matrix for the implicit scheme
mainDiagonal = 2 * (1 + alpha) * ones(numGridPoints, 1);
offDiagonal = -alpha * ones(numGridPoints, 1);
A = spdiags([mainDiagonal, offDiagonal, offDiagonal], [0 -1 1], numGridPoints, numGridPoints);
identityMatrix = speye(numGridPoints);
A([1 numGridPoints], :) = identityMatrix([1 numGridPoints], :); % No-flux boundary condition
% Define initial conditions and plot
sigma = domainLength / 16;
u = (1 / (sigma * sqrt(2 * pi))) * exp(-0.5 * (x / sigma).^2);
u = u';
plot(x, u, 'r'); hold on;
xlabel('$x$', 'Interpreter', 'latex', 'FontSize', 14);
ylabel('$u(x, t)$', 'Interpreter', 'latex', 'FontSize', 14);
title('Solution of the Diffusion Equation', 'Interpreter', 'latex', 'FontSize', 16);
% solve the linear system and plot at appropriate intervals
for timeStep = 1:numTimeSteps
b = [0; ...
alpha * u(1:numGridPoints-2) + 2 * (1 - alpha) * u(2:numGridPoints-1) + alpha * u(3:numGridPoints); ...
0];
u = A \ b;
if mod(timeStep, plotFrequency) == 0
plot(x, u, 'b');
end
end
Output:
This code effectively simulates the diffusion of a substance in a 1D domain, capturing the spreading of the initial Gaussian profile over time. The Crank-Nicolson method ensures stability and accuracy, making it suitable for long time simulations.
For more documentation on the functions used, refer the following links:
- ‘speye’: https://www.mathworks.com/help/matlab/ref/speye.html
- ‘spdiags’: https://www.mathworks.com/help/matlab/ref/spdiags.html
Hope this helps!
