NUMERICAL SOLUTIONS OF DIFFRENTIAL EQUATIONS
USING FAST FOURIER TRANSFORM
OLGA MARTIN
Department of Mathematics
University “Politehnica” of Bucharest
Splaiul Independentei 313,Bucharest 16
ROMANIA
Abstract. The numerical algorithm presented in this paper uses the Fast Fourier Transform (FFT) to find the coefficients of a trigonometric interpolation polynomial. Thus, a considerable reduction of the computations is achieved.
We illustrate this algorithm by a numerical example. First, we solve a Cauchy problem for an elastic vibrating system, using the finite difference method. Then, with the values of the approximate solution obtained in the equidistant points from the interval [0, 1], we shall find an interpolation polynomial using FFT. Also, we study the approximation of the numerical solution and stability of the difference scheme, which correspond of a second-order differential equation.
Key – Words: Fast Fourier Transform, elastic vibrating system, interpolation polynomial, approximation, stability.
1
1 Introduction
Fast Fourier Transform is a direct method for the solution a difference system. Many authors,[ 1],
[5], [9], [11] tackle now this modified Fourier method, which relies on the J.W. Cooley and J.W.Tucky algorithm.
In this paper on presented two numerical examples, which emphasize the construction of the interpolation polynomial of the solution of a difference problem and a numerical method to obtain the periodical solution of a Cauchy problem for a partial differential equation.
2 Problem formulation
Let on the interval [0, 2], 2N equidistant points, x0 = 0, x1,..., x2N-1, x2N= 2. A periodical function f :RR, with the period T = 2, attain the given values, fj = f (xj), at the points xj,
j0, 1,...., 2N - 1.
The interpolation trigonometric polynomial of the function f is defined by the equality
(1)
where
(2)
To halve the computations number for calculation coefficients Ak, Bk on separates the components of fj with even index from those with odd index and on defines
yk = f2k + i f2k+1, k0,1,...., N - 1 (3)
With Cooley and Tukey algorithm [11], we calculate the sums of the form
,
j0,1,...,N-1 (4)
where
(5)
The relations between Cjand Aj, Bj are given by
the next theorem, [11].
Theorem
(6)
(7)
where j0,1,...., N, B0 = BN, CN = C0.
The relations (4) can be written in a matriceal form
Y = TC (8)
where
and for ,
(9)
is a square matrix. Since: on
obtains:, where IN is the identity matrix
of order NxN and is the conjugate matrix of T. It
follows from the form of T that is a nonsingular matrix,
so
(10)
Thus, the relations (5) become
(11)
which define the Inverse Discrete Fourier Transform. Considering that we get
(12)
where P is a permutation matrix. From the relations (6)
we find Ai, Bi and the interpolation polynomial (1).
3Interpolation of reticular function
Consider a steel beam which is simply supported at the ends of length l = 25 cm and diameter d = 8mm. At his half is placed a weight G = 3 daN, which is acted by o force varying harmonically with time:
F(t) = F0.sin pt,
F0 = 5 daN and p = 100 rad/sec.
The elasticconstant
.
The equation of motion then becomes
(13)
where: g – the acceleration of gravity. The circular frequency is defined by
(14)
and the period is
. (15)
Multiplying (13) with ratio P/g we obtain, (16)
and introducing the numerical values we get Cauchy problem:
(17)
The complete solution of the problem (17) is composed of two functions.
The first of these represents a natural vibratory motion
(18)
The second vibratory motion is due to exciting force F(t)
(19)
with the period:. (20)
The interpolation polynomial will be the sum
P(t) = P1(t) + P2(t) (21)
where P1, P2 are the polynomials which correspond to u1, u2 respectively.
1. For the natural vibratory motion with T1 = 0,03 we change the interval 0, T1 in the interval 0, 1 by the relation: x = t/T1
.
Let on the interval [0, 1], N = 23 equidistant points with the step of division h = 1/8.
We solve the problem:
(22)
by finite difference method and the results of calculations are entered in the appropriate rows of Table 1:
Table 1
1
i / 0 / 1 / 2 / 3 / 4 / 5 / 6 / 7xi / 0 / 1/8 / 2/8 / 3/8 / 4/8 / 5/8 / 6/8 / 7/8
ui / 0 / -0.017 / -0.024 / -0.017 / 0 / 0.017 / 0.024 / 0.017
1
According with (3), we introduce the notation
yk = u1, 2k + i.u1, 2k+1 , k = 0,1, 2, 3.(23)
and because: , we define the permutation matrix .(24)
From equation (12)
where Y =
Therefore
(25)
and from (7)
A0 - iB0 = 0; A1 - iB1 = -0.0002 + 0.02445 i ;
A2 - iB2 = 0; A3 - iB3 = 0,0002 - 0,00005 i.
The interpolation polynomial of u1(x) has the form:
or in variable t
(26)
2. In the second vibratory motion with the period T2 = 0.06 sec the interval 0, 1 is substituted for 0, T2 by the relation: x = t/T2.
Using the same division for 0, 1 with the step h = 1/8, as in preceding case, we obtain for Cauchy problem:
(27)
the next values for u(xi), i0,1,…,7 by finite difference method:
Table 2
1
i / 0 / 1 / 2 / 3 / 4 / 5 / 6 / 7xi / 0 / 1/8 / 2/8 / 3/8 / 4/8 / 5/8 / 6/8 / 7/8
ui / 0 / 0.038 / 0.05 / 0.038 / 0 / -0.038 / -0.05 / -0.038
1
Then
yk = u2, 2k + i.u2, 2k+1 , k = 0,1, 2,
and
A0 - iB0 = 0; A1 - iB1 = -0.052 i ; A2 - iB2 = 0;
A3 - iB3 = - 0.0019 i.
The interpolation polynomial of the function u2(x) is
and in variable t
(28)
We obtain finally
4 Approximation of solution
In accordance with the difference scheme for the equations (22) and (27) we have:
and from Taylor’s formula
hence
where and C4(0,T) is the set of the function with the derivatives to four order continuous.
From Darboux property of function u there is 1, 2, for which
hence
Thus, the finite difference method used here is accurate to the order of h2.
5 Stability
In order to demonstrate this property for the difference problem which corresponding, for example, to (27), we consider in the interval [0, 1], N = 23 equidistant points, tn, n = 0, 1, 2,…, N, t0 = 0, tN = 1. Let be this partition of [0, 1] with the step h and we shall denote: un = u(tn), B = 2 and n = qsin ptn..
The difference problem may be written in the form:
(29)
where
(30)
In accordance with [3],[4] the scheme in difference (29) is stable, if for any f there is a unique solution and
(31)
where C is an independent constant of h.
An equivalent schema for (29) is:
(32)
Now we define
(33)
Let Y be the bi-dimensional space with un,
nY and we define the norm:
(34)
Since , we shall define a norm, which depends on h-step of partition of the interval [0, 1]:
(35)
Relationship between these norms is:
(36)
where: .
We shall prove that
(37)
for any linear operator T. Indeed,
Then, in accordance with [3] we obtain from the definition of the norm in Yh and (36):
.
and thus, (37) is true. Also, we have:
(38)
Recall that for any matrix T, [ 3]:
we have
(39)
Hence, if
or
and from (38):
because Nh = 1. From the formula (33)
We define
and
Then,
(40)
because: .
Therefore, there is a constant C = 2.eB such that (31) is fulfilled and the difference scheme for Cauchy problem (27) is stable.
References:
[1] BRIGHAM E.O., The Fast Fourier Transform,
Prentice-Hall, Englewood Cliffs, 1974.
[2] BHASKAR A., Criticality of damping in multi-degree-
of –freedom systems, J.of Apll. Mech., 64,
pp. 387-393,(1997) .
[3] GODUNOV S. K., REABENKI V.S., The theory of
difference schemes, Editura Tehnica, Bucharest,
1977.
[4] MARCIUK G. I., SHAYDUOV V., Refining of the
solutions of the differene schemes, Mir Publishers,
Moscow, 1983.
[5] MARCIUK G. I, Numerical Methods, Romanian Academia Publishing House, Bucharest, 1983.
[6] MARTIN O., Une methode de resolution de l’equation
dutransfert des neutrons, Rev. Roum. Sci. Tech.
Mec. Apll., 37, (1992), 623 - 649.
[7] MARTIN O., OLARIU V., Variational method for the
calculation of solutions of the neutron transfer equation,
Politehn. Univ.Bucharest Sci.Bull., Ser.D Mech. Engrg,
55, (1993), 5-12.
[8] PACHPATTE B.G., On some finite difference inequalities
in two independent variable, J. Math. Anal. Appl., 254,
(2001), 587-598.
[9] ROMBALDI J. E., Problemes corriges d’analyse
numerique, Masson, Paris, 1996.
[10] STAREK L., INMAM D., J., A symmetric inverse
vibration problem for non-proportional under-damped
system, J. of Appl. Mech., 64, (1997),601 – 605.
[11] STOER J., Introduction in numerical mathematics,
Springer – Verlag, Berlin, 1972.
1