According to the 12-step program, it’s time to solve a 1D diffusion equation, aka Heat equation of the form,

\displaystyle \frac{\partial u}{\partial t} = \nu \frac{\partial^2 u}{\partial x^2},

where {\nu} is the viscosity. The process is very similar to solving a 1D convection equation, but this time we need to approximate a second derivative. An easy way to get an expression that approximates a second derivative is by using Taylor series:

\displaystyle u(x+h,t)=u(x,t)+hu'(x,t)+\frac{h^2}{2!}u''(x,t)+\cdots,

and

\displaystyle u(x-h,t)=u(x,t)-hu'(x,t)+\frac{h^2}{2!}u''(x,t)+\cdots.

Now add them up and truncate to get,

\displaystyle u(x+h,t)+u(x-h,t)=2u(x,t)+h^2u''(x,t).

Since we are intereseted in the second derivative, we rewrite the last expression as

\displaystyle u''(x,t)=\frac{1}{h^2}\left[u(x-h,t)-2u(x,t)+u(x+h,t)\right].

The only part of the code that needs to change is the finite-difference loop itself:

// Finite-difference loop:
  for (int it=1; it<=nt-1; it++)
    {
      for (int k=0; k<=nx-1; k++)
	{
	  un[k][it-1] = u[k][it-1];
	}
	  // Keep values of u at the 'ends' constant:
      u[0][it] = un[0][it-1];
      u[nx-1][it] = un[nx-1][it-1];
      for (int i=1; i<=nx-2; i++)
	{
	  u[i][it] = un[i][it-1] + vis*dt/dx/dx*(un[i+1][it-1]-2*un[i][it-1]+un[i-1][it-1]);
	}
    }

Imagine that we have a rod of length 2m where the initial temperature distribution is the step-function as defined in the code. If no heat is added, we expect the rod to cool over time. The shown plot supports this.

 

About these ads