Computational Physics - Department of Physics

10.2 Diffusion equation 313

vec u(n+1); // This is uinAu = r

vec r(n+1); // Right side of matrix equation Au=r

// setting up the matrix

a = c = - alpha;

b = 2 + 2*alpha;

// Time iteration

for(int t = 1; t <= tsteps; t++){

// Calculate rforuse intridag, right hand side of the Crank Nicolson method

for(int i = 1; i < n; i++){

r(i) = alpha*u(i-1) + (2 - 2*alpha)*u(i) + alpha*u(i+1);

}

r(0) = 0;

r(n) = 0;

//Thensolve the tridiagonal matrix

tridag(a, b, c, r, u, xsteps+1);

u(0) = 0;

u(n) = 0;

// Eventualprintstatements etc

}

10.2.4Solution for the One-dimensional Diffusion Equation

It cannot be repeated enough, it is always useful to find caseswhere one can compare the
numerical results and the developed algorithms and codes with closed-form solutions. The
above case is also particularly simple. We have the following partial differential equation

a(t)

t

g(x)

b(t)

x

ui− 1 ,j+ 1 ui,j+ 1 ui+ 1 ,j+ 1

ui− 1 ,j ui,j ui+ 1 ,j

✲

✻

Fig. 10.3Calculational molecule for the Crank-Nicolson scheme.