Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

51, 2, pp. 251-264, Warsaw 2013

TREFFTZ FUNCTIONS FOR NON-STATIONARY PROBLEMS

Krzysztof Grysa, Beata Maciejewska

Kielce University of Technology, Department of Management and Computer Modeling, Kielce, Poland

e-mail: grysa@tu.kielce.pl

Different types ofTrefftz functions fornon-stationary linear andweaklynonlineardifferential
equations are presented. The Trefftz methods are defined and briefly described. Certain
results for non-stationaryproblemsof heat conduction (amongothers boundary temperature
identification and thermal diffusivity estimation), for beam vibration, for thermoelasticity
and for the wave equation (direct and inverse problem of membrane vibrations) are shown.
Inmany cases, the FEMwithTrefftz functions (FEMT) as probe functions is applied. Three
kinds of FEMT are tested for direct and inverse non-stationary problems. Examples of the
making use ofT-functions for solving inverse problems and smoothing inaccurate input data
are discussed.

Key words: non-stationary inverse problems, Trefftz method, FEMT

1. Introduction

The method known as “Trefftz method” was firstly presented in 1926 (Trefftz, 1926). An ap-
proximate solution of a problem is assumed to have form of a linear combination of functions
that satisfy the governing differential linear equation (without sources). In the case of using
finite elements, such an approximation is applied to each element. The unknown coefficients are
determined based on approximate fulfilling the boundary, initial and other conditions. It leads
to a system of algebraic equations. Such an approach to the approximate solving of a direct
or inverse problem is referred to as indirect Trefftz method (Kita and Kamiya, 1995). In order
to apply the Trefftz method, a complete system of the so-called T-functions (or T-complete
functions or Trefftz base) has to be known.
Generally, Trefftz bases fall into two broad classes. TheF-Trefftzmethod constructs its basis

function space by allowing many source points outside the domain of interest to keep the basis
functions regular. For each source point, the fundamental solution is adopted; therefore, this
method is also called themethodof fundamental solutions (MFS) (Liu et al., 2006).TheT-Trefftz
bases consist of functionswhich fulfill identically the considereddifferential equation,(Ciałkowski
and Frąckowiak, 2000; Grysa, 2010; Kita and Kamiya, 1995; Zieliński, 1995). Here we consider
only the T-Trefftz bases.
Starting at the end of the 70s, the T-functions have been presented for both bounded and

unbounded regions, for Laplace and Helmholtz equations (Herrera, 1984; Herrera and Sabina,
1978), biharmonic equation (Gurgeon and Herrera, 1981) and Stokes problem (Herrera and
Gurgeon, 1982). For elasticity problems,T-funtionshavebeen shownbyJiroušek andTeodorescu
(1982), Qin (2000). The basic sets of T-functions for stationary problems can be found by
Zieliński (1995).
The non-stationary problem was simplified by getting rid of the time variable. A result of

such an approach was the Helmholtz type equation.
In the late 90s, the T-functions for non-stationary problems have been published, at first,

based on the paperRosenbloom andWidder (1956), for parabolic equations (Futakiewicz, 1999;



252 K. Grysa, B. Maciejewska

Grysa, 2003; Hożejowski, 1999; Rakoczy, 1997), then for wave equation (Ciałkowski and Frąc-
kowiak, 2000; Ciałkowski and Jarosławski, 2003; Maciąg, 2007), beam vibration (Al-Khatib et
al., 2008), and for plate vibration and thermoelasticity (Maciąg, 2009b).
TheT-functions for the linear differential equationusuallyhave beenobtainedby themethod

of variable separation. In 2000, two new methods of deriving the T-functions were published
(Ciałkowski and Frąckowiak, 2000). The first one based on the expanding on the function into
Taylor series is particularly simple and effective; the second is based on the so-called inverse
operators (Ciałkowski and Frąckowiak, 2000, 2003).
During the last 10 years, theTrefftzmethod has been frequently applied to findapproximate

solutions of the direct and inverse problemsGrysa, 2010).
At first, a “global” approach was applied with good results for simple geometry and initial-

boundary conditions. Next, the FEM with Trefftz functions as trial functions was investigated.
Recently, the T-functions have been tested to smooth inaccurate internal data for the wave and
for heat transfer equations with very good results (Grysa and Leśniewka, 2010; Maciąg, 2009a).
Here we present the basic sets of T-functions for chosen nonstationary linear and weakly

nonlinear problems in bounded regions. Then, the Trefftz methods are briefly described with
a special attention to the least square method in relation to nonstationary direct and inverse
problems.Moreover, approximate solutions of problems of heat conduction andwavemotion are
shown.

2. Basic sets of T-functions for chosen differential equations

We confine our attention to non-stationary linear differential equations that model physical
phenomena in bounded areas. Sets of T-functions for stationary problems are discussed e.g. in
Ciałkowski andGrysa (2010), Zieliński (1995), and in themonograph byGrysa (2010). The sets
of T-functions presented below are useful in approximate solving of direct and inverse problems.
Below, in almost all cases (except for thermoelasticity problems), the governing equation are
presented in dimensionless form.

2.1. Heat conduction equation

In 1D, for (x,t) ∈ Ω× (0,T), Ω ⊂ R1, the T-functions (also called the heat polynomials)
read (Hożejowski, 1999; Rosenbloom andWidder, 1956)

vn(x,t) =
[n/2]∑

k=0

xn−2ktk

(n−2k)!k!
(2.1)

where [x] denotes the function floor(x).
In the polar system of coordinates (r,t), the T-complete set of functions consists of two sets

of functions

gn(r,t)=
n∑

k=0

(
r
2

)2n−2k
tk

((n−k)!)2k!
un(r,t)=−gn(r,t)lnr+ qn(r,t) (2.2)

with

qn(r,t) =
n∑

k=0

an−k
(
r
2

)2n−2k
tk

((n−k)!)2k!



Trefftz functions for non-stationary problems 253

and an defined by formulas

a0 =0 an =1+
1
2
+ . . .+

1
n
for n> 0

In 2D, T-functions are defined as (Ciałkowski and Frąckowiak, 2000)

Vm(x,y,t)= vn−k(x,t)vk(y,t) n=0,1, . . . k=0, . . . ,n (2.3)

and m=n(n+1)/2+k; vn(u,t) is given by (2.1).
In the cylindrical system of coordinates (r,z,t), the T-complete set of functions consists of

two sets of functions

V1m(r,z,t) = gn−k(r,t)vk(z,t) V2m(r,z,t) =un−k(r,t)vk(z,t) (2.4)

with n=0,1, . . .; k=0, . . . ,n and m=n(n+1)/2+k.
In 3D, we have (Ciałkowski and Frąckowiak, 2000)

Vnkl(x,y,z,t) = vn−k−l(x,t)vk−l(y,t)vl(z,t) (2.5)

Here n=0,1, . . .; k=0, . . . ,n; l=0, . . . ,k and vn(u,t) is given by (2.1).

2.2. Wave equation

In 2D, for (x,y,t) ∈ Ω × (0,T), Ω ⊂ R2, the variable separation method leads to the
following recurrent formulas for the T-functions (wave polynomials) Pmn0(x,y,t), Pmn1(x,y,t),
Qmn0(x,y,t) and Qmn1(x,y,t) (Maciąg, 2007)

P000 =1 Q000 =0

P(n−k)k0 =−
1
n

(
xQ(n−k−1)k0+yQ(n−k)(k−1)0+ tQ(n−k−2)k1+ tQ(n−k)(k−2)1

)

P(n−k−1)k1 =−
1
n

(
xQ(n−k−2)k1+yQ(n−k−1)(k−1)1+ tQ(n−k−1)k0

)

Q(n−k)k0 =
1
n

(
xP(n−k−1)k0+yP(n−k)(k−1)0+ tP(n−k−2)k1+ tP(n−k)(k−2)1

)

Q(n−k−1)k1 =
1
n

(
xP(n−k−2)k1+yP(n−k−1)(k−1)1+ tP(n−k−1)k0

)

(2.6)

for n=1,2, . . ., k=0,1, . . . ,n. If any subscript is negative, theT-functionwith such a subscript
is equal to zero.
In 3D, for n = 1,2, . . ., k = 0,1, . . . ,n, l = 0,1, . . . ,n and m = 0,1 and for P0000 = 1,

Q0000 = 0 the T-functions P(n−k−l−m)klm(x,y,z,t) and Q(n−k−l−m)klm(x,y,z,t) may be obta-
ined from the following recurrent formulas (Maciąg, 2009b)

P(n−k−l)kl0 =−
1
n

(
xQ(n−k−l−1)kl0+yQ(n−k−l)(k−1)l0+zQ(n−k−l)k(l−1)0

+ tQ(n−k−l−2)kl1+ tQ(n−k−l)(k−2)l1+ tQ(n−k−l)k(l−2)1
)

P(n−k−l−1)kl1 =−
1
n

(
xQ(n−k−l−2)kl1+yQ(n−k−l−1)(k−1)l1+yQ(n−k−l−1)k(l−1)1

+ tQ(n−k−l−1)kl0
)

Q(n−k−l)kl0 =
1
n

(
xP(n−k−l−1)kl0+yP(n−k−l)(k−1)l0+zP(n−k−l)k(l−1)0

+ tP(n−k−l−2)kl1+ tP(n−k−l)(k−2)l1+ tP(n−k−l)k(l−2)1
)

Q(n−k−l−1)kl1 =
1
n

(
xP(n−k−l−2)kl1+yP(n−k−l−1)(k−1)l1+yP(n−k−l−1)k(l−1)1

+ tP(n−k−l−1)kl0
)

(2.7)



254 K. Grysa, B. Maciejewska

2.3. Beam vibration equation

For the beam vibration equation

( ∂4

∂x4
+
∂2

∂t2

)
u=0 (x,t)∈Ω× (0,T) Ω⊂R1

the T-function set can be derived considering the Taylor series for a solution of the equation. In
order to eliminate the derivatives with respect to variable x of the order higher than four, the
following form of the governing equation will be used

∂4νu

∂x4ν
=(−1)ν

∂2νu

∂t2ν
ν =1,2, . . .

The sums of coefficients accompanying the same derivatives stand for the T-functions. Fi-
nally, the T-functions for the beam vibration equation for n = 1,2, . . . read (Al-Khatib et al.,
2008)

Sn(x,t)=
[An/2]∑

j=0

(−1)jtAn−2jxBn+4j

(An−2j)!(Bn+4j)!
(2.8)

where

An =n−3
[n−2
4

]
−3 Bn =4

[n−2
4

]
−n+5

2.4. Plate vibration equation

For the beam vibration equation

( ∂4

∂x4
+2

∂4

∂x2∂y2
+
∂4

∂y4
+
∂2

∂t2

)
u=0 (x,y,t)∈Ω× (0,T) Ω⊂R2

the T-functions Pkl(x,y,t) and Qkl(x,y,t) for n= 1,2, . . ., and for P00 = 1, Q00 = t may be
obtained from the following recurrent formulas (Maciąg, 2009b)

Pn1 = yP(n−1)0 Qn1 = yQ(n−1)0
Pn(n−1) =xP(n−1)(n−1) Qn(n−1) =xQ(n−1)(n−1)

P(n+2)0 =
2x(n+1)P(n+1)0−x

2Pn0−4t2P(n−1)0−2tQ(n−2)0
(n+1)(n+2)

Qn0 =
2x(n+1)Q(n−1)0−x

2Q(n−2)0−4t
2Q(n−4)0+2tPn0

(n+1)(n+2)

Pn(k+1) =
(n−k+1)Pn(k−1)−xP(n−1)(k−1)+yP(n−1)k

k+1
1¬ k¬n−1

Qn(k+1) =
(n−k+1)Qn(k−1)−xQ(n−1)(k−1)+yQ(n−1)k

k+1
1¬ k¬n−1

(2.9)

2.5. Thermoelasticity

The set of thermoelasticity equations read (Grysa andMaciąg, 2011; Maciąg, 2009b)

µ∇2u+(λ+µ)grad divu= ρü+γgradT
1
κ

∂T

∂t
=∇2T (2.10)



Trefftz functions for non-stationary problems 255

Here u represents the displacement vector, ∇ stands for nabla operator, µ, λ denote Lamé
constants, ρ means mass density, κ is the coefficient of thermal diffusivity. Applying Lamé
substitution u = gradφ+ rotΨ in Eq. (2.10)1 leads to the following wave equations for the
potentials φ and Ψ

(
∇2−

1
c21

∂2

∂t2

)
φ=mT

(
∇2−

1
c22

∂2

∂t2

)
Ψ=0 (2.11)

with c21 =(λ+2µ)/ρ, c
2
2 =µ/ρ standing for the velocities of waves and m= γ/(c

2
1ρ). After rew-

riting equations (2.11) in dimensionless form, theT-functions for themhave form (2.6) in 2Dand
form (2.7) in 3D. However, the T-functions are defined only for homogeneous equations. Hence,
for Eq. (2.11)1 one has to find a particular solution for known T . For heat transfer equation
(2.10)2, the T-functions are defined by (2.1), (2.3) or (2.5) in 1D, 2D and 3D, respectively.
Sets ofT-functions for stationaryproblemsarediscussed e.g. inCiałkowski andGrysa (2010),

Zieliński (1995) and in the monograph byGrysa (2010).

2.6. Weakly nonlinear problems

Consider a nonlinear partial differential equation Lu= f in Ω with the condition Bu= g
on a set of points K being a boundary ∂Ω in the case of the direct problem or a part of the
boundary and some points inside Ω in the case of an inverse problem. Let L=Llin+N; then
we have Llinu= f−Nuwith Llin standing for the linear part of the operator L and N being
its nonlinear part. We assume that the nonlinearity is weak what means that its influence on
the solution is, say, of the order of 5-10% of the linearized problem solution.
Assume that for the homogeneous equation Llinϕ = 0 the T-complete functions {ϕ

∗
n} are

known. Then the equation Llinu= f−Nu is a quasilinear equation in implicit form. To solve
the equation Llinu= f−Nu, the Picard iterations can be implemented

Llinu
(k) = f−Nu(k−1) in Ω (for k=0 Llinu

(0) = f) (2.12)

and

Bu(k) = g on K k=0,1, . . . ,K (2.13)

with u(k) = u
(k)
part+u

(k)
gen. Here u(k) is a linear combination of T-functions, k =0,1, . . . ,K. Its

form is obtained in the following way. Because u
(k)
gen is a general solution to the homogeneous

equation Llinu
(k)
gen = 0, it can be expressed as a linear combination of a finite number of the

T-functions {ϕ∗n}. Let us approximate the source function f and the nonlinear term Nu
(k−1)

as f ≈
∑
j bjPj and Nu

(k−1) ≈
∑
iciPi with Pi,Pj beingmonomials (f and Nu

(k−1) are then
approximated by polynomials). A particular solution to equation (2.12) is then

u
(k)
part =L

−1
lin

(∑

j

bjPj +
∑

i

ciPi
)

(2.14)

The stopping criterion for the Picard iteration can be adopted in the form

E = max
xm∈Ω

√√√√
M∑

m=1

(
u(k)(xm)−u(k−1)(xm)

)2
<ε (2.15)

with {xm}Mm=1 standing for a set of trial points; M and ε are chosen arbitrarily.
As an example, consider a 1Ddimensionless heat conduction problem in Ωwith anunknown

conductivity coefficient that depends on temperature



256 K. Grysa, B. Maciejewska

( ∂
∂x
λ(u)

∂

∂x
−
∂

∂t

)
u= f λ(u) =λ0+αu

α being a small number.
Assume that u

∣∣
∂Ω
is known, ∂Ω= {(x1, t),(x2, t)},Ω× (0, te),Ω⊂R. Then

L=
∂

∂x
λ(u)

∂

∂x
−
∂

∂t
=Llin+N

with

Llin =
∂2

∂x2
−
∂

∂t
N =

∂

∂x
αu
∂

∂x

For thehomogeneous equation, theT-functionsareknown, see (2.1).Ageneral solution to the
homogeneous equationhas thena form u(k)gen =

∑
na
(k)
n vn.With thepolynomial approximationof

the source term f andof thenonlinear term Nu(k−1), it is easy tofindaparticular solution u
(k)
part,

(2.14), to equation (2.12), making use of the following formulas

w0,k =L
−1
lin(t

k)=−
tk+1

k+1
w1,k =L

−1
lin(xt

k)=−
xtk+1

k+1

wm,k =L
−1
lin(x

mtk)=
1
k+1

(
m(m−1)wm−2k+1−x

mtk+1
)

k=0,1, . . . m=2,3, . . .

With u(k) = u
(k)
part + u

(k)
gen it is easy to find the coefficients a

(k)
n and – as a consequence

– a solution in the k-th iteration to equation (2.12) with conditions (2.13). Finally, stopping
criterion (2.15) shows whether the next iteration has to be calculated.
Convergence of the Picard iteration is here an open question. In the case when it fails, any

other iterating procedure can be adopted.

3. Trefftz methods

The Trefftz methods have not received a precise definition, although this terminology has had
wide acceptance. Herreras definition of what is meant by a Trefftz method is (Herrera, 1984):
Given a region of an Euclidean space or some partitions of that region, a “Trefftz Method”

is any procedure for solving initial boundary value problems of partial differential equations or

systems of such equations, on such region, using solutions of that differential equation or its

adjoint, defined in its subregions.
The Trefftz methods are generally divided into two groups: direct methods and indirect

methods. In the direct method, the governing equation is to be fulfilled in a weak formulation
with the trial function being a linear combination of T-functions. In the indirect method, the
solution itself is approximated as a linear combination of the T-functions.
Among the indirect methods the most popular ones seem to be the least square method

(LSM) and collocation method. Different formulation of the indirect methods are presented
in Kita and Kamiya (1995). The LSM seems also to be useful when dealing with the inverse
problems.
Let us consider a problem formulated as follows:

Lu= f in Ω× (0, te)

u= g1 on SD× (0, te) p
∂u

∂n
= g2 on SN × (0, te)

p
∂u

∂n
+αu= g3 on SR× (0, te)

u= g4 on Sint×Tint u=h on Ω for t=0



Trefftz functions for non-stationary problems 257

where Lu = f is a parabolic linear differential equation, Sint stands for a set of points inside
the considered region (SD∪SN ∪SR)⊂ ∂Ω, Tint ⊂ (0, te) is a set of moments of time and the
functions gi, i=1,2,3,4 and h are of proper class of differentiability in the domains in which
they are determined. In order to solve such a problem in an approximate way, let us consider
the following objective functional

I(v) =
∫

Ω×(0,te)

(Lv−f)2 dΩdt+w1
∫

SD×(0,te)

(v−g1)
2 dSdt+w2

∫

SN×(0,te)

(
p
∂v

∂n
−g2
)2
dSdt

+w3
∫

SR×(0,te)

(
p
∂v

∂n
+αv−g3

)2
dSdt+w4

∑

Sint×Tint

(v−g4)
2+w5

∫

Ω

(v−h)2 dΩ
(3.1)

Here wi, i=1, . . . ,5, stand for constant positive weighting parameters which are to preserve
the numerical equivalence between the terms, and v ∈ V . The problem may be formulated as
follows: find a function u∈V such that

I(u)= min
∀v∈V
I(v) (3.2)

Such an approach leads to the LSM.
Let V be a space generated by a finite number of T-functions for the equation Lu = 0. If

v is a linear combination of the Trefftz functions for the homogeneous equation plus vpart and
vpart is known (the function vpart stands for a particular solution to equation Lu= f), the first
term in equation (3.1) vanishes. If f =0 then the function vpart is equal to zero.
Usually, the weighting parameters are equal to one. Moreover, the objective functional may

be completed with a term or terms regularizing the solution (Ciałkowski et al., 2007; Grysa et
al., 2008).
The advantage of the LSM is that the stiffnessmatrix in the system of algebraic equations is

always symmetric and positive definite (or semi-definite) and that even complicated constraints
of the solution can be easily incorporated with a numerical method. However, in the case of
global approach, thedisadvantage ofLSMis that the conditionnumberwill increase significantly,
compared with FEMLi et al., 2007).
When the region Ω×(0, te) has to be divided into subregions, one introduces the time-space

elements (Grysa et al., 2009b; Maciejewska, 2004). Then the FEMwith T-functions (FEMT) is
convenient to be applied. At least three kinds of FEMT may be used:
(a) FEMT with the condition of continuity of temperature in the common nodes of elements

(b) FEMT with no temperature continuity at any point between the elements

(c) Nodeless FEMT (substructuring). Instead, in each finite element the temperature is ap-
proximated with a linear combination of the T-functions.

Hence, depending of the kind of FEMT, the T-functions are used as trial functions or the
approximate solutionbeinga linear combinationof theT-functions is constructed ineach element
separately (Grysa and Leśniewska, 2010). The elements may be greater than in the classic FEM
because theapproximate solution fulfill thegoverning equation identically.Moreover, thenumber
of base functions can be smaller than in the classic FEM.
Usually, on the borders between the elements, one considers continuity of the approximate

solution and its normal derivative. Here, we only requireminimization of the difference of values
v+ and v− and of values ∂v/∂n+ and ∂v/∂n−. Therefore, the objective functional is completed
with terms

∫

Γ×(0,te)

(v+−v−)2 dSdt
∫

Γ×(0,te)

( ∂v
∂n+
−
∂v

∂n−

)2
dSdt



258 K. Grysa, B. Maciejewska

accompanied with positive weights. Here Γ stands for the border between the elements v+

and v− (n+ and n−) denote values of the function v (normal vector n) on both sides of Γ ,
respectively. Equation (3.2) always leads to a system of algebraic equations for the coefficients
of the T-functions combination.
In the case of operator L being a hyperbolic one the analysis is similar.

4. Examples of the problems solved using Trefftz method

4.1. 2D heat conduction problem

Consider a transient 2D heat conduction problem in dimensionless coordinates in a square.
On three sides of the square, the heat flux is prescribed. The fourth thermal condition is known
inside the square on a line distant at δb ∈ (0,0.99) from the fourth side of the square. For
δb =0 the initial boundary problem is considered (a direct problem). For δb > 0 the considered
problem is called an inverse one. Hence, the heat conduction equation is considered

∆T =
∂T

∂t
(x,y)∈Ω=(0,1)× (0,1) t∈ (0, tK) (4.1)

with the conditions

T(x,y,t)
∣∣∣
t=0
=T0(x,y) T(x,y,t)

∣∣∣
x=0
=h1(y,t)= ey+2t

∂T

∂y
(x,y,t)

∣∣∣∣
y=1

=h2(x,t)= ex+1+2t
∂T

∂y
(x,y,t)

∣∣∣∣
y=0

=h3(x,t)= ex+2t

T(1− δb,yi, tk)=Tik

(4.2)

Here Tik stand for temperatures measured in points distant at δb ∈ (0,0.99) from the side
x=1 of the square, i=1, . . . ,I; k=1, . . . ,K.
In order to solve the problem, the time interval is divided into subintervals. In each subin-

terval, the domain Ω is divided into 4 finite subdomains Ωj, j = 1,2,3,4 and in each subdo-
main the temperature is approximated as a linear combination of the T-functions Vm(x,y,t),
Tj(x,y,t) ≈ T̃j(x,y,t) =

∑N
m=1 c

j
mVm(x,y,t), given by formula (2.3). In the case of the first

time subinterval, the initial condition is known. For the next subintervals, the initial condition
is understood as the temperature field in the subdomain Ωj at the final moment of time in
the previous subinterval. The LSM is used tominimize the inaccuracy of the approximate solu-
tion on the boundary, in the initial moment of time and on the borders between the elements.
This way the unknown coefficients of the combination are calculated. Moreover, on the border
between the elements the heat flux jumps are minimized (Grysa and Leśniewska, 2010).
The objective functional reads

J =
∑

i

∫

Di

(
(T̃i(x,y,0)−T0(x,y)

)2
dD+

∑

i

tm∫

tm−1

∫

Γi

(
T̃i(0,y,t)−h1(y,t)

)2
dΓ dt

+
∑

i

tm∫

tm−1

∫

Γi

(∂T̃i
∂y
(x,1, t)−h2(x,t)

)2
dΓ dt+

∑

i

tm∫

tm−1

∫

Γi

(∂T̃i
∂y
(x,0, t)−h3(x,t)

)2
dΓ dt

+
∑

i,j

tm∫

tm−1

∫

Γij

(T̃i− T̃j)
2 dΓ dt+

I∑

i

K∑

k=1

(
T̃j(1− δb,yi, tk)−Tik

)2

+
∑

i,j

tm∫

tm−1

∫

Γij

(∂T̃i
∂x
−
∂T̃j
∂x

)2
dΓ dt+

∑

i,j

tm∫

tm−1

∫

Γij

(∂T̃i
∂y
−
∂T̃j
∂y

)2
dΓ dt

(4.3)



Trefftz functions for non-stationary problems 259

where Γij is a common boundary of the elements Ωi and Ωj. The dimensionless time in-
terval is (0,0.01], i.e. t0 = 0. The temperatures are known at points (1 − δb,yi) with
yi ∈ {0.1,0.2,0.3,0.4,0.6,0.7,0.8,0.9} for t ∈ {0.0025,0.0050,0.0075,0.01} and δb ∈ (0,0.99).
It means that I =8, K =4 and m=1,2,3,4. For δb =0, the problem becomes a direct one.
Temperatures inside the square (the internal temperature responses, abbr. ITRs) were genera-
ted from the accurate solution and then interfered with a random noise N(0,0.05). The noisy
“measurements” have been next approximated with the linear combination of the T-functions
due to the formula

T(1− δb,y,t)≈ T̃δb(y,t)=
S∑

s=1

asVs(1− δb,y,t) (4.4)

For the calculation, we have taken S=18T-functions. Because the exact solution is known
(Grysa and Leśniewska, 2010), the error of the approximate solution in the square in the i-th
time subinterval can be described using the norm

δL2 =

√√√√√√√

∫

D

(
T̃(x,y,t)−T(x,y,t)

)2
dD

∫

D

(
T(x,y,t)

)2
dD

·100% D=Ω× (ti−1, ti) (4.5)

The error on the part of the boundary where the temperature has been calculated in the
m-th time subinterval reads

δL2
∣∣∣
x=1
=

√√√√√√√√√

tm∫
tm−1

dt
1∫
0

(
T̃(1,y,t)−T(1,y,t)

)2
dy

tm∫
tm−1

dt
1∫
0

(
T(1,y,t)

)2
dy

·100% (4.6)

The results obtained with the nodeless FEMT for noisy and smoothed data are presented in
Table 1.

Table 1.Norm δL2 and δL2
∣∣
x=1
as functions of δb and ITRs

Noisy ITRs Smoothed ITRs Exact ITRs
δb δL2 [%] δL2|x=1 [%] δL2 [%] δL2|x=1 [%] δL2 [%] δL2|x=1 [%]

0.1 0.7204 1.4844 0.04318 0.0745 0.0315 0.07938
0.3 1.0455 0.6235 0.04470 0.0833 0.0478 0.08107
0.5 0.6929 0.4371 0.04069 0.0991 0.0469 0.07197
0.7 1.1008 0.9198 0.04463 0.0932 0.0459 0.08005
0.9 0.7263 0.4261 0.04339 0.0890 0.0464 0.07696
0.99 0.4265 0.3864 0.04573 0.0884 0.0517 0.08119

Comparing the errors of the unknown boundary identification for the exact and smoothed
noisy data, one observes a slight difference between their values. For the exact ITRs, the results
are sometimes a bitworse than those for the smoothed data. Itmay be a result of approximating
the solution with a small number of T-functions.

4.2. Thermal diffusivity estimation

In Grysa et al. (2009a), thermal diffusivity estimation based on temperature measurements
in PTFE (teflon), ice and soil are reported. In the experiments, a thermal probe having a shape



260 K. Grysa, B. Maciejewska

of a thin tube with 16 sensing or heating elements inside was used. The probe was placed in a
specimen of shape of a hollow cylinder with outer radius much greater than the outer radius of
the probe. Detailed technical specifications referring to the probe can be found in Spohn et al.
(2007). In each specimen every 310s, one of the sensors was used as a heater (heating interval
lasted 900s) and next the responses were registered with 10s sampling (Grysa et al., 2009a).
The coefficients of thermal diffusivity ofPTFE, ice and soil were known, and the experimentwas
to verify the method of thermal diffusivity estimation in the specimen with unknown thermal
properties (Sendek, 2007).
When the probe is in themode of heating, the heat is dissipated to the surroundingmedium

in the radial direction. The dimensionless heat transfer equation in cylindrical coordinates

∂2T

∂ξ2
+
1
ξ

∂T

∂ξ
+
∂2T

∂ς2
=
∂T

∂τ

ξ∈ (ξp,ξb) ς ∈ (0, ςk) τ > 0
(4.7)

with the boundary condition T(ξb, ς,τm) = Tm is considered. Here ξp and ξb stand for the
dimensionless outer radius of theprobe (theborderbetween theprobeand the specimen)and the
outer radius of the specimen, respectively, ςk denotes length of the k-th sensor, τ =κt/(rp−rb)
is the dimensionless time.Moreover, κ denotes the thermal diffusivity, rp and rb stand for outer
radius of the probe and outer radius of the specimen, respectively.
The problem of thermal diffusivity estimation was solved using T-functions. For cylindrical

geometry, the T-functions for the heat conduction equation are presented in Section 2, see
formula (2.2). An approximate solution θ(ξ,ς,tau) to the problem formulated above (note that
only one boundary condition is prescribed and no initial condition is formulated!) has form of
a linear combination of T-functions. The unknown coefficients of the combination are chosen to
minimize the functional

J =

ςk∫

0

M∑

m=1

(
θ(ξp, ς,τm)−Tm

)
dς (4.8)

for k=1,2, . . . ,16. Here M denotes the number of temperature measurements. Mimization of
the functional J means that we demand θ(ξp, ς,τ) to fit the boundary condition as accurately
as possible. The necessary condition for minimum of the functional leads to a system of alge-
braic linear equations for coefficients of the linear combination of T-functions, standing for the
approximate solution θ(ξ,ς,τ).
The functional J can be calculated only for a given value of thermal diffusivity (here denoted

as κ). In further consideration, we focus our attention on how J varies with change of κ. Our
objective is to find such κ for which J(κ) is minimal. First we need to choose a range of
the parameter κ bracketing the minimum. Next, in order to find such avalue of κ for which
the functional J has the minimum value, a method called the golden section is used (Press
et al., 2007). The process of searching the value of κ continues until the bracketing interval is
acceptably small. Of course, the smaller the sampling step for κ, the better the accuracy of
the thermal diffusivity. Finally, the value of thermal diffusivity for the investigated medium (a
sample) is taken as a mean value for the results obtained for each sensor.
In the considered problem, only 12T-functions have been taken into account. Three intervals

of heating were considered. The results are presented in Table 2.
As it can be seen, the obtained diffusivities are at least acceptable when compared to the

reference values, presented in the first column. It is worth noting that such results have been
achieved for only 12 T-functions and for only one boundary condition taken into account.



Trefftz functions for non-stationary problems 261

Table 2.Average values of thermal diffusivity

Medium / κ value 1st interval 2nd interval 3rd interval

PTFE / κ=1.13 ·10−7 1.23 ·10−7 1.89 ·10−7 1.62 ·10−7

Ice / κ=1.6 ·10−6 nomeasurements 1.58 ·10−6 1.77 ·10−6

Soil / κ=8.72 ·10−7 4.77 ·10−7 6.30 ·10−7 6.26 ·10−7

4.3. Free vibrations of a square membrane

Consider free vibrations of a squaremembrane Ω=(0,1)×(0,1) described by dimensionless
wave equation in 2D (Grysa, 2010; Maciąg, 2009a), with the conditions

u(x,y,0)=x(x−1)y(y−1)
∂u

∂t

∣∣∣∣
(x,y,0)

=0

u(0,y,t) =u(x,0, t) =u(x,1, t) = 0

(4.9)

and with the displacement measured in three inner points (the internal responses, abbr. IRs)

u(1−ε,0.25, t) =u1(t) u(1−ε,0.50, t) =u2(t) u(1−ε,0.75, t) =u3(t) (4.10)

If ε> 0wehave an inverse problem (themeasureddisplacements are then called the internal
responses, abbr. IRs) andwe search for the solution in thewhole domainbut especially u(1,y,t).
In the case of ε=0, we consider a direct (initial-boundary) problem.
The exact solution to the problem for u(1,y,t) = 0 is known (Maciąg, 2009a), and is used to

generate the IRs. The solution is approximated in the whole square according to formula (4.4),
i.e. u≈ w =

∑N
n=1cnVn with Vn being T-functions (wave polynomials) for the wave equation

in 2D (presented in Section 2, formulas (2.6)). The objective functional has the following form

I =
1∫

0

1∫

0

(
w(x,y,0)−x(x−1)y(y−1)

)2
dxdy+

1∫

0

1∫

0

(∂w
∂t
(x,y,0)

)2
dxdy

+
∆t∫

0

1∫

0

(
w(0,y,t)

)2
dydt+

∆t∫

0

1∫

0

(
w(x,0, t)

)2
dxdt+

∆t∫

0

1∫

0

(
w(x,1, t)

)2
dxdt

+
∆t∫

0

{ 3∑

k=1

(
w(1−ε,0.25k,t)−uk(t)

)2}
dt

(4.11)

The approximation of the solution with the linear combination of wave polynomials of the
order up to 7 shows that the difference between the exact solution for the vibration as a function
of time for the location x=1, y=0.5 and the approximation for ε=0, ε=0.05 and ε=0.3
does not exceed 0.2% in the case when IRs are undisturbed.
Now let us disturb the internal response ukl =u(1−ε,0.1(k−1),0.1(l−1)), k,l=1,2, . . . ,11

according to the formula udkl =ukl(1+ξkl)where ξkl are randomnumbers of normal distribution
with the expected value equal to zero and standard deviation equal tp 0.01 (N(0,0.01)). The
disturbance ξkl has been generated separately for each ukl using a random-number generator.
In this case, the relative error

E =

√√√√√√√√

∆t∫
0

(
u(1,0.5, t)−w(1,0.5, t)

)2
dt

∆t∫
0

(
u(1,0.5, t)

)2
dt

(4.12)

has a value from circa 2000% for ε=0 (direct problem) to circa 180000% for ε=0.8.



262 K. Grysa, B. Maciejewska

Itmeans that avery small disturbance causes averybadapproximation.However, smoothing
thedisturbeddatabya linear combination ofwavepolynomials in a similarway like in example 1
(formula (4.4)) leads to much better results. The coefficients of the combination are chosen so

as tominimize
∑
kl

(
wε(1−ε,(k−1)0.1,(l−1)0.1)−udkl

)2
. After smoothing, the error is similar

as for the exact data (Maciąg, 2009a).

5. Final conclusions

The object of this paper is to present the existing T-function bases for bounded regions for non-
stationary problems.Theapplications ofT-functions for solving non-stationary inverse problems
show the usefulness of the functions. Many others are reported by Grysa (2010). The use of T-
functions to solve weakly nonlinear direct and inverse problems is a new idea. A conclusion
that follows the presented examples (and that should be checked bymathematicians) is that an
approximate solution of a direct or inverse problem with the use of T-functions seems to be of
better quality if the functions describing conditions and input data are formulated in the same
subspace of the space generated by the Trefftz base. Increasing the number of T-functions leads
initially to better results (Grysa and Leśniewska, 2010), but usually their number cannot be
too large, because it leads to ill-conditioned problems. Therefore, the best seems to use them as
basic functions in the FEM, since their number in an element is usually not too large.
There is still a lot of open problems left here. They concern some technical problems as well

as mathematical background of the Trefftz methods.
For instance, when dealing with the FEM or FEMT with time-space elements, a question

concerning the dependencebetween the length of the spatial and the “length” of the time dimen-
sions of the element appears. Probably theymay depend on the investigated signal propagation
velocity. In the case of the heat conduction problems, the “velocity” of temperature propagation
is difficult to define. Theoretically, it is equal to infinity, but practically, it depends strongly on
the accuracy of temperature measurements.
Generally, in FEMT, big finite elements can be used. However, their size depends on the

number of trial functions. The kind of the dependence is not known.
In the second example1 (and in many papers devoted to the inverse heat conduction pro-

blems), small numbers of points with input data (with ITRs) lead to good results. Also the
incomplete data (e.g. lack of initial conditions (Sendek, 2007), or boundary conditions known
only on a part of the boundary (Maciąg, 2009b)) lead to good (comparable with the accurate)
results.
Many other questions arise when applying Trefftz methods. From the point of view of engi-

neers, obtaining better results with an increasing number of T-functions is a good justification
to apply the method. However, the mathematicians are kindly invited to check some ideas of
the engineers. It would be useful to unify the studies of the two communities.

References

1. Al-Khatib M.J., GrysaK., Maciąg A., 2008, Themethod of solving polynomials in the beam
vibration problems, Journal of Theoretical and Appled Mechanics, 46, 347-366

2. Ciałkowski M.J., Frąckowiak A., 2000,Heat Functions and Their Applications to solving the
Heat Conduction and Mechanical Problems, Wyd. Politechniki Poznańskiej [in Polish]

3. Ciałkowski M.J., Frąckowiak A., 2003, The heat functions and related in solving chosen
problems ofmechanics,Part I – Solving somedifferential equationswith the use of inverse operator,
Studia i Materiały, Technika, 3, 7-70, Uniwersytet Zielonogórski



Trefftz functions for non-stationary problems 263

4. Ciałkowski M.J., Frąckowiak A., Grysa K., 2007, Physical regularization for inverse pro-
blems of stationary heat conduction, Journal of Inverse and Ill-Posed Problems, 15, 347-364

5. Ciałkowski M.J., Grysa K., 2010, Trefftz method in solving the inverse problems, Journal of
Inverse and Ill-Posed Problems, 18, 595-616

6. Ciałkowski M.J., Jarosławski M., 2003, The use of symbolic computation for generating of
wave equations solution, Zeszyty Naukowe Politechniki Poznańskiej, M. R. i T., 56, 115-140 [in
Polish]

7. Futakiewicz S., 1999, Method of heat functions to solving direct and inverse problems of heat
conduction, PHDTheses, Politechnika Poznańska [in Polish]

8. Grysa K., 2003,Heat polynomials and their applications, Archives of Thermodynamics, 24,
107-124

9. Grysa K., 2010,Trefftz Functions and Their Application in Solving Inverse Problems, Wyd. Po-
litechniki Świętokrzyskiej [in Polish]

10. Grysa K., Hożejowski L., Marczewski W., Sendek-Matysiak E., 2009, Thermal diffusi-
vity estimation from temperature measurements with a use of a thermal probe, Proc. Int. Conf.
Experimental Fluid Mechanics, 63-71

11. Grysa K., Leśniewska R., 2010, Different finite element approaches for inverse heat conduction
problems, Inverse Problems in Science and Engineering, 18, 3-17

12. Grysa K., Leśniewska R., Maciąg A., 2008, Energetic approach to direct and inverse he-
at conduction problems with Trefftz functions used in FEM, Comp. Ass. Mech. Eng. Sci., 15,
171-182

13. Grysa K., Maciąg A., 2011, Solving direct and inverse thermoelasticity problems by means of
Trefftz base functions for finite element method, Journal of Thermal Stresses, 34, 1-16

14. Grysa K., Maciąg A., Maciejewska B., 2009,Wave polynomials as base functions of FEM in
problems of elastokinetics,Proc. of the 7th EUROMECH Solid Mech. Conf., 191-192

15. Gurgeon H., Herrera I., 1981, Boundary methods. C-complete systems for biharmonic equ-
ations, [In:] Boundary Element Methods, Brebbia C.A. (Edit.), New York, CML Publ. Springer,
431-441

16. Herrera I., 1984,BoundaryMethods. AnAlgebraic Theory, Boston,PitmanAdvancedPublishing
Program

17. Herrera I., 2000, Trefftz method: A general theory, Numerical Methods for Partial Differential
Equations, 16, 561-580

18. Herrera I., Gurgeon H., 1982, Boundary methods, c-complete systems for Stokes problem,
Computer Methods in Applied Mechanics and Engineering 30, 25-41

19. Herrera I., Sabina F., 1978, Connectivity as an alternative to boundary integral equations,
Construction of Bases Appl. Math. Phys. Sci., 75, 2059-2063

20. Hożejowski L., 1999, Heat polynomials and their application in direct and inverse problems od
heat conduction, PHDTheses, Politechnika Świętokrzyska [in Polish]

21. Jiroušek J., Teodorescu P., 1982, Large finite elements method for the solution of problems
in the theory of elasticity,Computers and Structures, 15, 575-587

22. Kita E., Kamiya N., 1995, Trefftz method: an overview,Advances in Engineering Software, 24,
3-12

23. Li Z.C., LuT.T., HuangH.T., ChengA.H.-D., 2007,Trefftz, collocation, and other boundary
methods – a comparison,Numerical Methods for Partial Differential Equations, 23, 93-144

24. LiuR.-F.,YeihW.,KuoS.-R.,ChenY.-W., 2006, IndirectT-Trefftz andF-Trefftzmethods for
solving boundary value problem of Poisson equation, Journal of the Chinese Institute of Industrial
Engineers, 29, 989-1006



264 K. Grysa, B. Maciejewska

25. Maciąg A., 2007,Wave polynomials in elasticity problems,Eng. Trans., 55, 129-153

26. Maciąg A., 2009a, The usage of wave polynomials in solving direct and inverse problems for
two-dimensional wave equation, International Journal for Numerical Methods in Biomedical Engi-
neering, DOI: 10.1002/cnm.1338

27. MaciągA., 2009b,Trefftz Functions for Chosen Direct and Inverse Problems ofMechanics,Wyd.
Politechniki Świętokrzyskiej [in Polish]

28. Maciejewska B., 2004, Application of the modified method of finite elements for identification
of temperature of a body heated with a moving heat source, Journal of Theoretical and Appled
Mechanics, 42, 771-788

29. PressW.H.,TeukolskyS.A.,VetterlingW.T., FlanneryB.P., 2007,Numerical Recipies,
The Art of Scientific Computing, Third Ed. Cambridge University Press

30. Qin Q.-H., 2000,The Trefftz Finite and Boundary Element Method, Boston: WIT Press, South-
ampton

31. Rakoczy K., 1997,On some class of the parabolic polynomials,Fasciculi Mathematici, 27, 81-93

32. Rosenbloom PC, Widder DV., 1956, Expansion in terms of heat polynomials and associated
functions,Transactions of the American Mathematical Society, 92, 220-266

33. Sendek E., 2007, Investigation of thermal properties of a medium with the use of thermal probe
Mupus-Pen, PHDTheses, Politechnika Świętokrzyska [in Polish]

34. SpohnT., SeiferlinK., HagermannA.,Knollenberg J., BallA.J., BanaszkiewiczM.,
Benkhoff Jo., Gadomski S., Gregorczyk W., Grygorczuk J., Hłond M., Kargl G.,

Kuhrt E., Komlee N., Krasowski J., Marczewski W., Zarnecki J.C., 2007, Mupus – a
thermal andmechanical properties probe for the rosetta lander philae,Space ScienceReviews,128,
339-362

35. Trefftz E., 1926,EinGegenstueck zumRitz’schenVerfahren,Proc. 2nd Int. Congress of Applied
Mechanics, 131-137

36. Zieliński A.P., 1995, On trial functions applied in the generalized Trefftz method, Advances in
Engineering Software, 24, 147-155

Funkcje Trefftza dla problemów niestacjonarnych

Streszczenie

W pracy przedstawione są różne rodzaje funkcji Trefftza dla niestacjonarnych liniowych i słabo nie-
liniowych równań różniczkowych. Zdefiniowane są i krótko opisane metody Trefftza. Prezentowane są
niektóre wyniki dla niestacjonarnych problemów przewodzenia ciepła (m.in. identyfikacja temperatury
na brzegu oraz estymacja dyfuzyjności termicznej), drgań belki oraz dla termosprężystości i równania
falowego (prosty i odwrotny problem drgańmembrany).W pracy pokazany jest również przykład zasto-
sowania metody elementów skończonych z funkcjami Trefftza (MEST) jako funkcjami próbnymi. Trzy
rodzaje MEST są testowane na prostych i odwrotnych problemach niestacjonarnych.W pracy omawia-
ne są również przykłady zastosowania T-funkcji w rozwiązywaniu problemów odwrotnych w połączeniu
z wygładzaniem niedokładnych danych wejściowych.

Manuscript received February 21, 2012; accepted for print May 14, 2012