One-dimensional heat equation

Let a one-dimensional heat equation with homogenous Dirichlet boundary conditions and zero initial conditions be subject to spatially and temporally distributed forcing

The second derivative operator with Dirichlet boundary conditions is self-adjoint with a complete set of orthonormal eigenfunctions, , . This information can be used to diagonalize operator which facilitates straightforward determination of the frequency response. For systems with spatially varying coefficients and non-normal generators the frequency response analysis is typically done numerically using finite dimensional approximations of the differential operators. For example, the pseudospectral method with collocation points can be used to transform the frequency response operator into an matrix. However, for systems with differential operators of high order, spectral differentiation matrices may be poorly conditioned and implementation of boundary conditions may be challenging. In our method, numerical approximation of differential operators in the evolution equation is avoided by first recasting the system into corresponding two point boundary value problems and then using state-of-the-art techniques for solving the resulting boundary value problems with accuracy comparable to machine precision.

Two point boundary value representations of , , and

Application of the temporal Fourier transform yields the two point boundary value representation of the frequency response operator ,

where

	—	imaginary unit
	—	temporal frequency
	—	second-order differential operator in
	—	point evaluation functional at the boundary .

The two point boundary value representation for the adjoint of the frequency response operator, , is given by

Fig. : Block diagram of .

The representation of the operator is obtained by combining the two point boundary value representations of and . As illustrated in Fig. , this can be achieved by setting in ,

Integral form of a differential equation

We next illustrate the procedure for converting the differential equation into its corresponding integral form. By introducing an auxiliary variable

and by integrating we obtain

Here, and denote the indefinite integration operators of degrees one and two, the vector contains the constants of integration which are to be determined from the boundary conditions, and

The integral form of the 1D heat equation is obtained by substituting into ,

Now, by observing that

we can use to express the constants of integration in terms of ,

Finally, substitution of into yields an equation for ,

This equation only contains indefinite integration operators and point evaluation functionals which are known to be well-conditioned.

Matlab codes

% The system parameters:
om = 0; % set temporal frequency to the value of interest

dom = domain(-1,1);     % domain of your function
y = chebfun('y',dom);
fone = chebfun(1,dom);      % one function
fzero = chebfun(0,dom);     % zero function

% The system operators can be constructed as follows:
A0{1} = [-1i*om*fone, fzero, fone];
B0{1} = -fone;
C0{1} = fone;

% The boundary condition matrices are given by:
Wa0{1} = [1, 0];   % 1*phi(-1) + 0*D^{(1)}*phi(-1) = 0
Wb0{1} = [1, 0];   % 1*phi(+1) + 0*D^{(1)}*phi(+1) = 0

% The singular values and the associated singular functions of the
% frequency response operator can be computed using the following code

Nsigs = 25;    % find the largest 25 singular values
LRfuns = 1;    % and associated left singular functions

[Sfun, Sval] = svdfr(A0,B0,C0,Wa0,Wb0,LRfuns,Nsigs);

% Analytical expressions for the singular values and corresponding singular
% functions are given by:
Sa = (4./(([1:Nsigs].*pi).^2)).';    % analytical singular values
Sf1 = sin((1/sqrt(Sa(1))).*(y+1));   % analytical solution of the first singular function
Sf2 = sin((1/sqrt(Sa(2))).*(y+1));   % analytical solution of the second singular function

% The absolute error of the first 25 singular values
norm(Sval - Sa)

ans =

   5.4288e-15

plot(1:Nsigs,Sval,'bx','LineWidth',1.25,'MarkerSize',10)
hold on
plot(1:Nsigs,Sa,'ro','LineWidth',1.25,'MarkerSize',10);
xlab = xlabel('n', 'interpreter', 'tex');
set(xlab, 'FontName', 'cmmi10', 'FontSize', 20);
h = get(gcf,'CurrentAxes');
set(h, 'FontName', 'cmr10', 'FontSize', 20, 'xscale', 'lin', 'yscale', 'lin');

Twenty five largest singular values of the frequency response operator.

hold off; plot(y,-Sfun(:,1),'bx-','LineWidth',1.25,'MarkerSize',10)
hold on; plot(y,Sf1,'ro-','LineWidth',1.25,'MarkerSize',10);
xlab = xlabel('y', 'interpreter', 'tex');
set(xlab, 'FontName', 'cmmi10', 'FontSize', 20);
h = get(gcf,'CurrentAxes');
set(h, 'FontName', 'cmr10', 'FontSize', 20, 'xscale', 'lin', 'yscale', 'lin');
axis tight

The principal singular function of the frequency response operator.

hold off; plot(y,Sfun(:,2),'bx-','LineWidth',1.25,'MarkerSize',10)
hold on; plot(y,Sf2,'ro-','LineWidth',1.25,'MarkerSize',10);
xlab = xlabel('y', 'interpreter', 'tex');
set(xlab, 'FontName', 'cmmi10', 'FontSize', 20);
h = get(gcf,'CurrentAxes');
set(h, 'FontName', 'cmr10', 'FontSize', 20, 'xscale', 'lin', 'yscale', 'lin');
axis tight

The singular function of the frequency response operator corresponding to the second largest singular value.