The typical example of a parabolic eqn is the Heat Equation,
No matter how `rough' the initial data for the heat equation, the solution is differentiable infinitely often. Because of the rapid decay of transients, the heat eqn is often used as a means of reaching a steady-state (i.e., the solution to the Laplace eqn).
Three principle methods for numerical finite differences:
Stability requires 
This requires the solution of a linear system of equations. Analysis shows unconditional stability.
This again requires the solution of a linear system of equations. Analysis shows unconditional stability.
Here is another perspective on the computation. Instead of using finite differences to approximate differentiation, lets mimic the physics. We believe the diffusion process is the result of particles experiencing many random `kicks', by other particles. These repeated kicks drive down the concentration, removing heat in the process. The code below uses a random walk to mimic the `kicks' - at each instant in time, each particle takes a step of size chosen from a Gaussian distribution with zero mean and specified variance. The variance translates into the diffusion constant. In this code, particles are reflected perfectly at the walls.
Run this code, and see how it compares with a finite difference solution. How does changing the variance effect the solution?
program heat2
c test code to solve heat eqn by random walk inside a tube
c PURE WALK. INTERPOLATION ONTO GRID ONLY AT END
c c==0 outside tube
c DESIGN TUBULE TO SIT BETWEEN x=0.2 and x=0.8 i.e. i=21 - i=81
implicit real*8(a-h,o-z)
dimension u(101),up(101),uex(101),upp(101),
1 conc(10000),pos(10000)
common/grid/x(101)
common/Igrid/nx
idum=-1957
width=4.0
npart=2000
nx=101
dx=1.0d0/(nx-1)
do 01 i=1,nx
x(i)=(i-1)*dx
01 continue
t=0.0d0
c viscosity
dnu=0.10d0
c DATA u is concentration.
do 10 i=1,nx
u(i)=0.0d0
up(i)=0.0d0
10 continue
u(51)=1.0d0/dx
uex(51)=u(51)
do 14 k=1,npart
pos(k)=0.50d0
conc(k)=(1.0d0/npart)*(1.0/dx)
14 continue
c PRINT
sum=0.0d0
do 16 i=1,nx
sum=sum+u(i)
16 continue
write(6,665)t,sum
665 format('','AT TIME T= ',f12.7,' ||u|| = ',f12.7)
c time stepping
dt=0.5*dx*dx/dnu
c want variance = 2*dnu*dt
var=(2.0d0*dnu*dt)
do while (t .lt. 0.20)
t=t+dt
c reflecting boundary conditions for walker
do 90 k=1,npart
step=pos(k)+gauss(idum,var)
if (step .le. 0.0d0) then
diff=0.d0-step
pos(k)=0.d0+diff
else if (step .ge. 1.d0) then
diff=step-1.d0
pos(k)=1.d0-diff
else
pos(k)=step
endif
90 continue
enddo
sum=0.0
do 93 k=1,npart
sum=sum+pos(k)
93 continue
expect=sum/npart
write(6,94) expect
94 format('','EXPECTATION POSITION ',f12.7)
sum=0.0
do 95 i=1,npart
sum=sum+(pos(i)-expect)**2
95 continue
VAR=sum/(npart)
write(6,96)VAR
96 format('','VARIANCE = ',f12.7)
do 100 i=1,nx
do 100 k=1,npart
up(i)=up(i)+conc(k)*d(width,pos(k),x(i))*dx
100 continue
c PRINT
sum=0.0d0
do 116 i=1,nx
sum=sum+up(i)
116 continue
write(6,665)t,sum
c EXACT SOLUTION EXPLICIT
s=0.0d0
ds=0.45*dx*dx
dsdx2=ds/(dx*dx)
do while (s .lt. t)
s=s+ds
do 122 i=2,nx-1
upp(i)=uex(i)+dnu*dsdx2*(uex(i+1)-2.0d0*uex(i)+uex(i-1))
122 continue
upp(nx)=uex(nx)+dnu*dsdx2*(uex(nx-1)-uex(nx))
upp(1)=uex(1)+dnu*dsdx2*(uex(2)-uex(1))
do 128 i=1,nx
uex(i)=upp(i)
128 continue
enddo
dmass=0.0d0
do 130 i=1,nx
dmass=dmass+uex(i)
130 continue
write(6,131) dmass
131 format('','MASS OF EXACT SOLUTION IS ',f12.7)
open(unit=1,file='out')
do 1111 i=1,nx
write(1,1166)x(i),up(i)
1111 continue
write(1,1167)
write(1,1168)
do 1112 i=1,nx
write(1,1166)x(i),uex(i)
1112 continue
write(1,1169)
1166 format('',2(f12.7,1x))
1167 format('','"upart"')
1168 format('',' ')
1169 format('','"uex"')
stop
end
C**************************************************************************
C FUNCTION SUBROUTINES
C**************************************************************************
function d(dlen,w,z)
implicit real*8(a-h,o-z)
common/Igrid/nx
dx=1.0d0/(nx-1)
pi=4.0d0*datan(1.0d0)
dist=dabs(w-z)
if (dist .gt. 0.50d0*dlen*dx) then
d=0.0d0
else
d=(1.0d0/(dlen*dx))
1 *(1.0d0+dcos(pi*dist/(0.50d0*dlen*dx)))
endif
return
end
C----------------------------------------------------------------------------
c RANDOM CHOICE FROM WITHIN A GAUSSIAN DISTRIBUTION
c WITH MEAN 0 AND VARIANCE sigma
function gauss(idum,sigma)
implicit real*8(a-h,o-z)
data iset /0/
if (iset .eq. 0) then
02 v1=2.0d0*rand(idum)-1
v2=2.0d0*rand(idum)-1
r2=v1*v1 + v2*v2
if (r2 .ge. 1.0d0) go to 02
fact=dsqrt(-2.0d0*(sigma)*dlog(r2)/r2)
gset=v1*fact
gauss=v2*fact
iset=1
else
gauss=gset
iset=0
endif
return
end
C----------------------------------------------------------------------------