## Parabolic equations

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:

1. EXPLICIT

Stability requires

2. IMPLICIT

This requires the solution of a linear system of equations. Analysis shows unconditional stability.

3. CRANK-NICHOLSON The previous two schemes are 1st order accurate in time. This is a 2nd order scheme.

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----------------------------------------------------------------------------```

E. Bruce Pitman
Wed Apr 7 10:53:29 EDT 1999