Onedimensional heat equation
Let a onedimensional 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
selfadjoint 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 nonnormal 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
stateoftheart 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 
 —  secondorder 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
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 ,
Note that svdfr.m utilizes the integral form of a two point boundary value representation of
. This procedure yields accurate results even for
systems with high order differential operators and poorlyscaled coefficients.
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 wellconditioned.
Matlab codes
om = 0;
dom = domain(1,1);
y = chebfun('y',dom);
fone = chebfun(1,dom);
fzero = chebfun(0,dom);
A0{1} = [1i*om*fone, fzero, fone];
B0{1} = fone;
C0{1} = fone;
Wa0{1} = [1, 0];
Wb0{1} = [1, 0];
Nsigs = 25;
LRfuns = 1;
[Sfun, Sval] = svdfr(A0,B0,C0,Wa0,Wb0,LRfuns,Nsigs);
Sa = (4./(([1:Nsigs].*pi).^2)).';
Sf1 = sin((1/sqrt(Sa(1))).*(y+1));
Sf2 = sin((1/sqrt(Sa(2))).*(y+1));
norm(Sval  Sa)
ans =
5.4288e15
Singular values of the frequency response operator: svdfr versus analytical results.
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.

The principal singular function: svdfr versus analytical results.
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.

The singular function corresponding to the second largest singular value:
svdfr versus analytical results.
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.

