Sources for the heat equation and its adjoint

Above shows how a heat flow operator transforms the input q_t to the output q_t+1. What we really want to do is completely different. We want to define the input as the source function and we want to define the output as some arbitrary subset of q_t, q_t+1, q_t+2, etc. The adjoint we seek will carry us from the specified subset of q_t back to an image $\hat s_t$ of the sources s_t.

Defining $H=[1+\sigma delta_{xx}]$(where delta_xx is the second difference operator) reduces equation (3) to  

$$q_{t} \quad =\quad H_{t} q_{t-1} + s_{t}$$ (5)

where each q_t is a function of x. Recursive use of equation (5) over three time levels gives the matrix form  

$${\bold H}^{-1} \, \bold q \quad =\quad \left[ \begin{array}... ...\ s_1 \\ s_2 \\ s_3 \end{array} \right] \quad =\quad\bold s$$ (6)

A recursive solution begins at the top with q₀=s₀ and propagates downward. The idea for computing is first to load up q_t with the sources $q_t\leftarrow s_t$and then add in the heat-flow equation solution as you move along with

$$q_{t} \quad +=\quad H_{t} q_{t-1} \quad\quad\quad {\rm increasing}\ t$$ (7)

The adjoint of H is H'. We define the image of the sources $\tilde s_t$ by the recursion

$$\tilde s_{t-1} \quad =\quad H'_{t} \ \tilde s_{t} + q_{t-1} \quad\quad\quad {\rm decreasing}\ t$$ (8)

which can be summarized as  

$$({\bold H}^{-1})' \, \tilde {\bold s} \quad =\quad \left[ \... ...q_1 \\ q_2 \\ q_3 \end{array} \right] \ \quad =\quad\bold q$$ (9)

Equation (6) implies $\bold q=\bold H \bold s$and equation (9) implies $\tilde {\bold s} = ((\bold H^{-1})')^{-1}\bold q$.Basic mathematics tells us that $(\bold H^{-1})'=(\bold H')^{-1}$ which means that $((\bold H^{-1})')^{-1}=\bold H'$so equation (9) implies $\tilde {\bold s} = \bold H'\bold q$.Thus the operators defined by recursive solutions of equations (6) and (9) are adjoints.

As before, the computing methodology is to preload $\tilde s_{t}\leftarrow q_t$ and then use C-style augmenting assignments

$$\tilde s_{t-1} \quad +=\quad H'_{t} \ \tilde s_{t}$$ (10)

subroutine heatsteps( adj, sigma, ss, nx,nt,  qq )
integer               adj,            nx,nt,         it
real		           sigma, ss( nx,nt), qq( nx,nt)
temporary real	     	          tt( nx,nt)
call    adjnull( adj, 0,          ss, nx*nt,  qq, nx*nt)
if( adj == 0) {
	call copy( nx*nt, ss, tt)
	do it= 2, nt
		call heatstep( 0,1, sigma, tt(1,it-1), nx, tt(1,it))
	call copy( nx*nt, tt, qq)
} else {
	call copy( nx*nt, qq, tt)
	do it= nt,2,-1 
		call heatstep( 1,1, sigma, tt(1,it-1), nx, tt(1,it))
	call copy( nx*nt, tt, ss)
	}
return; end

The problem with the above subroutine is that it wastes memory, assuming that we can store the temperature field over all space and time. A more realistic assumption is that we can store it over all space but only at two time levels. Moving the copy commands into the inner loops, we require memory for only the field at time t and time t+1.

subroutine heatsteps( adj, sigma, ss, nx,nt,  qq )
integer               adj,            nx,nt,         it
real		           sigma, ss( nx,nt), qq( nx,nt)
temporary real			                        tt( nx)
call    adjnull( adj, 0,          ss, nx*nt,  qq, nx*nt)
if( adj == 0)  {
		call copy( nx, ss(1, 1), qq(1,1))	       # load initial
	do it= 2, nt {
		call copy( nx, ss(1,it),		   tt) # load source
		call heatstep( 0,1, sigma, qq(1,it-1), nx, tt)
		call copy( nx, tt,         qq(1,it))	       # extract 
		}
} else {
		call copy( nx,   qq(1,nt),         ss(1,nt))   # load initial
	do it= nt-1,1,-1 {
		call copy( nx,   qq(1,it), tt)		       # load 
		call heatstep( 1,1, sigma, tt, nx, ss(1,it+1))
		call copy( nx,   tt,               ss(1,it))   # extract source
		}
	}
return; end

Notice that the input ss(,) and output qq(,) still fill all of time and space so they appear to be impractically large. What saves us is that in practice there will not be sources everywhere in time and space nor will there be recording of the thermal field everywhere either. Thus in practice the copy calls will be reduced to copying the application operator inputs and outputs which should be much smaller.