How to solve diffusion equation by the crank - Nicolson method?

I have a diffusion equation 1D:

dC/dt =D*d2C/dx2

with D is changed with time.

Hepl me......

Thank you so much!

Answer 1

[Prasanna]

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!

[Torsten]

D is said to depend on t.