next up previous
Next: About this document Up: PDEs Previous: Hyperbolic equations

Parabolic equations

The typical example of a parabolic eqn is the Heat Equation,

displaymath781

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

    displaymath783

    Stability requires tex2html_wrap_inline785

  2. IMPLICIT

    displaymath787

    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.

    displaymath789

    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