The projective–iterative method and neural network estimation for buckling of elastic plates in nonlinear theory

Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088 www.elsevier.com/locate/cnsns

The projective–iterative method and neural network estimation for buckling of elastic plates in nonlinear theory Aliki D. Muradova

a,1

, Georgios E. Stavroulakis

b,c,*

a

b

Institute of Mechanics, Department of Mathematics, University of Ioannina, Greece Department of Production Engineering and Management, The Technical University of Crete, University Campus, GR-73100 Chania, Greece c Technical University of Braunschweig, Germany Received 3 June 2003; received in revised form 10 September 2005; accepted 10 September 2005 Available online 27 October 2005

Abstract The paper deals with a nonlinear buckling problem for von Karman elastic plates in bending. The simply supported and partially clamped plates are considered. A variational-projective approach with an iterative scheme is suggested for the calculation of eigenfunctions and the buckling loads for the studied problem. The convergence of the projective–iterative method is investigated. The bifurcation scenario is demonstrated by examples. Numerical results show the effectiveness of the proposed method. Prediction of the eigenvalues (which are the bifurcation points of the nonlinear problem) of the linearized problem with different sizes of the plate can also be done by an approximately trained neural network, as it is briefly demonstrated in this work.  2005 Elsevier B.V. All rights reserved. AMS classification: 74G60; 74G15; 74K20 Keywords: Elastic plates; Bifurcation; von Karman theory; Variational-projective method; Iterative scheme; Neural networks

1. Introduction In this paper the nonlinear problem of elastic rectangular plates in bending is considered. In the post-buckling area multiplicity of the solution and bifurcation points exist. This phenomenon for von Karman plate in bending has been investigated by many authors [1–5,8,9,12,13,20,21,23]. Here a numerical method for the calculations of buckling loads for rectangular plates is proposed. The method is based on a variational principle and uses iterations for the calculation of the eigenfunctions of the considered problem. As the variational approach we take the Galerkin projective method, where the *

Corresponding author. Address: Department of Production Engineering and Management, The Technical University of Crete, University Campus, GR-73100 Chania, Greece. Tel.: +30 28210 37418; fax: +30 28210 69410. E-mail address: [email protected] (G.E. Stavroulakis). 1 Present address: Department of Mineral Resources Engineering, The Technical University of Crete, Chania, Greece. 1007-5704/$ - see front matter  2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2005.09.001

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

1069

coordinate functions are chosen in accordance to the boundary conditions. The convergence of the projective – iterative scheme is investigated. The numerical results show branches of solutions and dependence of their number on the intensity of loading efforts k. In the last section we propose a neural network for the approximation of the buckling loads [22,24,27–29] (see also the related work on inverse problems reported in [25,26]). A back propagation neural network is applied to approximate the nonlinear relation between the plate dimensions l1, l2, and the corresponding buckling loads k = [k1, . . . , km]. The suggested algorithm can also be applied for more complicated plate buckling problems and for general estimation problems, arising in continuum mechanics [11,13,21]. 1.1. Formulation of the problem The system of equations which govern the von Karman plate bending and stretching problem reads: DD2 w  2h0 ½w; w þ kLw ¼ 0 E D2 w þ ½w; w ¼ 0; 2 2

2

2

2

ðx; yÞ 2 G; ð1Þ 2

2

w o w where ½w; w ¼ ooxw2 ooyw2 þ ooxw2 ooyw2  2 oxo oy , w denotes the deflection of the plate, w is the Airy stress function, ox oy which gives the stresses along the middle surface of the plate (see, e.g., [6]), w, w 2 W2,2(G), (W2,2(G) is the  ¼ ½0; l1   ½0; l2 , Sobolev 3 space of order (2, 2)), G is an open subset in R2, occupied by the plate, G 2Eh0 D ¼ 3ð1m2 Þ is the cylindrical rigidity of the plate, E is the Joungs modulus, h0 the width of the plate, m is the Poissons ratio, k denotes in physical sense an intensity of efforts (middle on the thickness of the plate) and L = D, o2/ox2. A detailed analysis for the existence of the solutions of (1) can be found, among others, in paper [3], monograph [7], etc. Following these results let

w ¼ 0;

L1 w ¼ 0;

w ¼ 0;

L2 w ¼ 0

ð2Þ

denote the boundary conditions on the edges of the plate. According to the known results if k 6 k1, where k1 is the first eigenvalue of the linearized problem DD2 w þ kLw ¼ 0

ðx; yÞ 2 G;

ð3Þ

w ¼ 0;

ðx; yÞ 2 oG;

ð4Þ

L1 w ¼ 0

then (1) and (2) has a unique trivial solution. If k > k1, then bifurcation takes place and (1) and (2) has at least three solutions (w, w), (w, w) (w 5 0, w 5 0) and (0, 0). We will study cases, when the plate is simply supported and partially clamped. When the plate is simply supported the boundary conditions read: ðaÞ ðbÞ

w ¼ 0; Dw ¼ 0 w ¼ 0; Dw ¼ 0

on oG;

ð5Þ

or for the stress function more appropriate condition in physical sense can be considered (see, e.g., Schaeffer and Golubitsky [20]) ow o ¼ ðDwÞ ¼ 0 on oG. on on If the plate is partially clamped then 8 < ow ¼ 0 for x ¼ 0; l1 ; w ¼ on : Dw ¼ 0 for y ¼ 0; l2 ;

ð6Þ

ð7Þ

w ¼ Dw ¼ 0 on oG; ð8Þ ow o ¼ Dw ¼ 0 on oG. ð9Þ on on Conditions (7), (8) and (7), (9) correspond to the clamped sides x = 0, l1 and simply supported ends y = 0, l2, respectively.

1070

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

For the calculation of the solutions in the case k > k1 the projective–iterative scheme is suggested. First the projective method is applied to (1) and (2). Next, after getting the system of nonlinear algebraic equations by the iterative method on the basis of the incremental loading method and Newton iterations a search of eigenfunctions is realized. First, we will give an analysis of existence for the bifurcation points for both cases L = D and L = o2/ox2. When L = o2/ox2 the bifurcation scenario for partially clamped plate has been studied, for instance, by Chien et al. [4]; Chien et al. [5]. In the most general case the bifurcation phenomenon was studied numerically by Allgower and Georg [1]; Caloz and Rappaz [8]; Crouzeix and Rappaz [9], etc. Besides, Holder and Schaeffer [12] and Schaeffer and Golubitsky [20] have shown that the so-called mode jumping, when the primary solution branches lose stability through further bifurcations, can occur only in the case of the partially clamped plate. Returning to Eqs. (1) for simplicity of mathematical analysis we take dimensionless coefficients in (1)–(3). 2. The bifurcation scenario for (1) and (2) 2.1. The case L = D For the plate with simply supported conditions the problem (3) and (4) has non trivial solutions, when k are eigenvalues of Laplacian, i.e., kmn = p2((m/l1)2 + (n/l2)2) and correspondingly these solutions are eigenfunctions of Laplacian. Thus, values kmn will be bifurcation points of solutions of the nonlinear problem (1) and (2). Obviously, some of them can be multiple bifurcation points. For instance, for simply supported square plate (l1 = l2) double bifurcation points are kmn = knm, m 5 n. When the plate is partially clamped we apply the procedure, described in the proceedings [4,5] (this procedure has been used for L = o2/ox2). For this purpose we find the solution of the linearized problem (3) and (7). According to the rule of separation of variables the solution is presented as wðx; yÞ ¼ uðxÞvðyÞ 6¼ 0. Substituting it into the governing equation (3) yields uðivÞ v þ 2u00 v00 þ vðivÞ u þ kðu00 v þ v00 uÞ ¼ 0.

ð10Þ

The boundary conditions (7) correspondingly will be u0 ð0Þ ¼ u0 ðl1 Þ ¼ 0;

uð0Þ ¼ uðl1 Þ ¼ 0;

vð0Þ ¼ vðl2 Þ ¼ 0;

v00 ð0Þ ¼ v00 ðl2 Þ ¼ 0.

y : n 2 N for the space of functions We choose the basis sin pn l2 fv 2 C 4 ½0; l2  : vð0Þ ¼ v00 ð0Þ ¼ vðl2 Þ ¼ v00 ðl2 Þg. Then (10) reduces into a fourth-order linear differential equation for u uðivÞ þ ðk  2ðpn=l2 Þ2 Þu00 þ ððpn=l2 Þ4  kðpn=l2 Þ2 Þu ¼ 0

ð11Þ

with the boundary conditions uð0Þ ¼ u0 ð0Þ ¼ uðl1 Þ ¼ u0 ðl1 Þ ¼ 0. Now, if we consider a solution of the form u = e 4

2

2

4

ð12Þ ax

then the characteristic equation reads

2

a þ ðk  2ðpn=l2 Þ Þa þ ððpn=l2 Þ  2 Þ Þ ¼ 0. qkðpn=l ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi The roots of this equation are ±pn/l2,  ðpn=l2 Þ2  k. For the studied problem only the case k > (pn/l2)2 has a physical sense. Therefore, the general form of the solution will be uðxÞ ¼ C 1 ea1 x þ C 2 ea1 x þ C 3 cos a2 x þ C 4 sin a2 x; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 where a1 = pn/l2, a2 ¼ k  ðpn=l2 Þ . Hence using the boundary conditions we obtain     1 a2 1 a2 C3 þ C4 ; C1 ¼  C4  C3 C2 ¼ 2 2 a1 a1

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

1071

and the system with respect to C3 and C4.   a2 ðcos a2 l1  cosh a1 l1 ÞC 3 þ sin a2 l1  sinh a1 l1 C 4 ¼ 0; a1  ða2 sin a2 l1 þ a1 sinh a1 l1 ÞC 3 þ a2 ðcos a2 l1  cosh a1 l1 ÞC 4 ¼ 0. This system has non trivial solution if the determinant of the coefficients equals zero, i.e.,       k k 2a2 þ ea1 l1 a1  sin a2 l1  a2 cos a2 l1  ea1 l1 a1  sin a2 l1 þ a2 cos a2 l1 ¼ 0. 2a1 2a1 Solving this equation with respect to k we define the bifurcation points (see Fig. 1 with n = 1). Table 1 demonstrates the dependence of k on the shape of the plate. In Table 1 the second eigenvalue k21, n = 1 of (3) and (7) is given. The first one for any n = 1, 2, . . . for different l1 with fixed l2 will be the same, for 2 2 example, when l2 = 1 it 0 is easily calculated by k1n ¼ pl2n . 2

2.2. The case L = o2/ox2



In this case for the simply supported plate, kmn ¼

p2 lm1

þ



nl1 ml2

2 2

and multiple bifurcation points depend

on l1, l2. If the conditions (7) are fulfilled then following the results of [4,5] the bifurcation points can be found from the solution of the equation (see Fig. 2) sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! pffiffiffi pffiffiffi 2p2 n2 k 4p2 n2 cos kl1  cos k  2 l1 ¼ 0. ð1  cos kl1 Þ þ 2 2 l2 l2

250 1 2

200

150

100

50

0

–50

–100

–150

–200 0

50

100

150

200

250

Parameter lambda

Fig. 1. Definition of the buckling loads: (1) solid line for l1 = 1, l2 = 1; (2) dashed line for l1 = 2, l2 = 1.

1072

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

Table 1 The eigenvalue k21 of (3) and (7) (L = D) for various plate dimensions l1, l2

1

2

3

4

5

1 2 3 4 5

37.7996 14.6174 11.5877 10.733 10.3854

37.748 9.4499 4.95873 3.65436 3.14163

38.54 9.29613 4.19996 2.59607 1.94012

38.9132 9.43699 4.12033 2.36248 1.61177

39.1051 9.55442 4.15131 2.31853 1.51198

5

6

x 10

4

2

0

–2

–4 1 2 –6

0

50

100

150

200

250

300

Parameter lambda

Fig. 2. Definition of the buckling loads: (1) solid line for l1 = 1, l2 = 1; (2) dashed line for l1 = 2, l2 = 1.

Table 2 The eigenvalue k21 of (3) and (7) (L = o2/ox2) for various plate dimensions l1, l2

1

2

3

4

5

1 2 3 4 5

66.5526 47.8394 43.49 41.8144 41.0008

44.8756 16.6382 13.2615 11.9599 11.1585

41.763 12.4278 7.39474 6.07729 5.84969

40.741 11.2189 5.8797 4.15954 3.5188

40.2798 10.7065 5.28253 3.44578 2.6621

pffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi kðkþ2mÞ 4p2 n2 ðkþmÞ2 4n2 p2 If we suppose k  l2 l1 ¼ 2mp, kl1 ¼ 2ðk þ mÞp, then k0 ¼ kðkþ2mÞl2 , where l1 ¼ will be a double n 2

2

bifurcation point. Table 2 gives the second eigenvalue (n = 1) of (3) and (7) for different sizes of the plate. The 2 2 first eigenvalue for n = 1, 2, . . . with different l1 and fixed l2 is the same. k1n ¼ 4pl2n , when l2 = 1. 2 For the approximation of the eigenvalues for other plate dimensions one can use neural network [22,24,27– 29] as it is demonstrated at the end of this paper.

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

1073

3. The projective method for (1) and (2) We find the solution of (1) and (2) by expanding into the double Fourier series ([20]): W N ðx; yÞ ¼

N X

wijN xij ðx; yÞ;

WN ðx; yÞ ¼

i;j¼1

N X

wijN uij ðx; yÞ;

ð13Þ

i;j¼1

where xij, uij are coordinate functions. Naturally, xij, uij are chosen in accordance to the boundary conditions. By applying the Galerkin method for (1) and (2) we have ðD2 W N ; xmn Þ ¼ ð½W N ; WN ; xmn Þ  kðLW N ; xmn Þ; ðD2 WN ; umn Þ ¼ ð½W N ; W N ; umn Þ;

m; n ¼ 1; 2; . . . ; N ;

where (Æ, Æ) denotes a scalar product in L2. Hence, we obtain a system of nonlinear algebraic equations with mn respect to wmn N , wN : mn mn K mn 1 wN ¼ A1;N ðwN ; wN Þ  kB wN ;

K mn 2 wN

¼

Amn 2;N ðwN ; wN Þ;

ð14Þ

m; n ¼ 1; 2; . . . ; N .

ð15Þ

Here K1, K2 are linear stiffness matrices (K1 is the bending and K2 is the stretching stiffness matrix, respectively), A1,N, A2,N are nonlinear stiffness matrices, arising from the second-order geometric nonlinearity terms, and B is a linear matrix, arising after discretization of the operator L. By using the double Fourier series expansion for the solution of (1) and (2) we assume: 1 1 X X wij xij ðx; yÞ; Wðx; yÞ ¼ wij uij ðx; yÞ. ð16Þ W ðx; yÞ ¼ i;j¼1

i;j¼1

Then the system (14) and (15) is, in fact, the truncated system with respect to the infinite one: mn mn K mn 1 w ¼ A1 ðw; wÞ  kB w;

K mn 2 w

¼

Amn 2 ðw; wÞ;

ð17Þ

m; n ¼ 1; 2; . . . ; 1

ð18Þ

(we use the notation, Ai,1 = Ai). 3.1. Simply supported plate For the simply supported plate the eigenfunctions of the Laplacian are considered as basic functions, i.e., xij ðx; yÞ ¼ p2ffiffiffiffiffiffi sin pi x sin pj y, uij(x, y) = xij(x, y) (see [16]). The use of the eigenfunctions of Laplacian as a spel1 l2 l1 l2

cial basis for getting a priori estimates is presented in details in Lions monograph [15]. The method is, in addition, is applied for some bending problems in [11]. As we have shown before the eigenvalues of (3) and (5) will    2  2 2 m be kmn ¼ p þ ln2 and these numbers are bifurcations points for the simply supported plate. In that l1    2 2 2   mn mn 4 m þ ln2 , A1,N = A2,N and their elements are case K 1 ¼ K 2 ¼ kmn , kmn ¼ p l1 Amn N ðwN ; wN Þ ¼

N X

wijN wkl N aði; k; m; j; l; nÞ;

i;j;k;l¼1

aði; k; m; j; l; nÞ ¼ a1ikm

p2 5=2

h i 2 2 ða1ikm a1jln þ a2ikm a2jln Þðil  jkÞ  ða1ikm a2jln þ a2ikm a1jln Þðil þ jkÞ .

2ðl1 l2 Þ cos pðm þ i þ kÞ  1 cos pðm  i  kÞ  1 þ ; ¼ mþiþk mik

Bmn ¼ kmn ;

when L ¼ D;

2

Bmn ¼ p2 ðm=l1 Þ ;

a2ikm ¼

cos pðm þ i  kÞ  1 cos pðm  i þ kÞ  1 þ ; mþik miþk

when L ¼ o2 =ox2 .

1074

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

If now for the stress function boundary conditions (6) are fulfilled, then uij ðx; yÞ ¼ p2ffiffiffiffiffiffi cos pi x cos pj y and l1 l2 l1 l2  mn mn again K 1 ¼ K 2 ¼ kmn , but Amn 1;N ðwN ; wN Þ ¼

N X

p4 ðl1 l2 Þ

5=2

wijN ðaði; j; m; nÞ  ijðm  iÞðn  jÞÞwNmi;nj

i;j¼1

þ aði; j; m; nÞðwNim;jn þ wNim;nj þ wmi;jn Þ  ijðm þ iÞðn þ jÞwNmþi;nþj N

ijðm  iÞðn þ jÞwNmi;nþj  ijðm þ iÞðn  jÞwNmþi;nj ; N h X p4 2 ij w ði2 ðj  nÞ  ijði  mÞðj  nÞÞwim;jn Amn N N 2;N ðwN ; wN Þ ¼ 5=2 ðl1 l2 Þ i;j¼1 þ ði2 ðj þ nÞ2  ijði  mÞðj þ nÞÞwim;jþn þ ði2 ðj  nÞ2  ijði þ mÞðj  nÞÞwiþm;jn N N i 2 iþm;jþn 2 þði ðj þ nÞ  ijði þ mÞðj þ nÞÞwN ; where a(i, j, m, n) = 0.5(i2(j  n)2 + j2(i  m)2). According to the obtained results, from (17) and (18) we get: mn mn kmn wmn N ¼ A1;N ðwN ; wN Þ þ kkmn w ;

ð19Þ

kmn wmn N

ð20Þ

¼

Amn 2;N ðwN ; wN Þ;

m; n ¼ 1; 2; . . . ; N

and, accordingly, the infinite system is mn kmn wmn ¼ Amn 1 ðw; wÞ þ kkmn w ;

ð21Þ

kmn wmn

ð22Þ

where kmn

Amn 2 ðw; wÞ;

¼ m; n ¼ 1; 2; . . . ; 1;   2 ¼ kmn ; pm . Let us estimate the coefficients wmn, wmn in (21) and (22). l1

Lemma 1. For the coefficients of the expansion of functions W(x, y) and W(x, y), which satisfy (1) and (5) or (1) and (5)(a), (6) for fixed k, the following estimations hold: wmn ¼ O



m 2 þ n2

2  ;

wmn ¼ O



m2 þ n2

2 

.

ð23Þ

Proof. Since, the series (16) converge, then wmn, wmn ! 0, when m, n ! 1. By virtue of (21) and (22), the mn definitions kmn , Amn 1 , A2 and kmn for fixed k we will have the formulae (23). h Let us rewrite (21) and (22) by the following way: ðkmn  kkmn Þwmn ¼ Amn 1 ðw; wÞ; wmn ¼ 

1 mn A ðw; wÞ; kmn 2

m; n ¼ 1; 2; . . .

ð24Þ ð25Þ

Hence, by virtue of the definitions of A1, A2 we obtain that the bifurcation points can be computed by solving the equation kmn  kkmn ¼ 0. 3.2. Partially clamped plate qffiffiffi In this case we use the following coordinate functions: xij ðx; yÞ ¼ X i ðx; l1 Þ l22 sin pj y, uij ðx; yÞ ¼ p2ffiffiffiffiffiffi l2 l1 l2 pj pj 2ffiffiffiffiffiffi pi p cos sin or u ðx; yÞ ¼ cos , where X (x, l ) satisfies (7) . As X (x, l ) we can take sin pi i 1 1 i 1 ij l l1 l2 l2 l1 l2 i qffiffiffih1 pði1Þ pðiþ1Þ pði1Þ pðiþ1Þ 2 cos x  cos x þ ði þ 1Þ sin x  ði  1Þ sin x . Then l1 l1 l1 l1 l1 m2;n mþ2;n K mn þ a2;mn wmn ; 1 wN ¼ a1;mn wN N þ a3;mn wN

ð26Þ

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

where 2

a1;mn ¼ ða2m1 þ b2n Þ ð1 þ ðm þ 1Þðm  3ÞÞ; a2;mn ¼ ða2m1 þ b2n Þ2 ð1 þ ðm þ 1Þ2 Þ þ ða2mþ1 þ b2n Þ2 ð1 þ ðm  1Þ2 Þ; 2

a3;mn ¼ ða2mþ1 þ b2n Þ ð1 þ ðm  1Þðm þ 3ÞÞ; pm pn am ¼ ; bn ¼ . l1 l2 N X Amn ðw ; w Þ ¼ wijN wkl N N N a1 ði; k; m; j; l; nÞ; 1;N i;j;k;l¼1

a1 ði; k; m; j; l; nÞ ¼

8p4 ðl1 l2 Þ7=2

ðaði; k; m; j; l; nÞ  2jklbði; k; m; j; l; nÞÞ;

for the boundary condition (8): n h i aði; k; m; j; l; nÞ ¼ l2 ði  1Þ2 I 1i1;k;m  ði þ 1Þ2 I 1iþ1;k;m þ ði þ 1Þði  1Þ2 I 2i1;k;m  ði  1Þði þ 1Þ2 I 2iþ1;k;m h io þðjkÞ2 I 1i1;k;m  I 1iþ1;k;m þ ði þ 1ÞI 2i1;k;m  ði  1ÞI 2iþ1;k;m ajln ; h i bði; k; m; j; l; nÞ ¼ ði  1ÞI 1k;i1;m þ ði þ 1ÞI 1k;i1;m þ ði þ 1Þði  1ÞI 3i1;k;m  ði  1Þði þ 1ÞI 3iþ1;k;m bjln ; I 1ikm ¼ bi;m1;k  bi;mþ1;k þ ðm þ 1Þdi;k;m1  ðm  1Þdi;k;mþ1 ; I 2ikm ¼ dm1;i;k  dmþ1;i;k þ ðm þ 1Þai;k;m1  ðm  1Þai;k;mþ1 ; I 3ikm ¼ di;k;m1  di;k;mþ1 þ ðm þ 1Þbi;k;m1  ðm  1Þbi;k;mþ1 ; when (9) is fulfilled: n h i 2 2 2 2 aði; k; m; j; l; nÞ ¼ l2 ði  1Þ I 3i1;k;m  ði þ 1Þ I 3iþ1;k;m þ ði þ 1Þði  1Þ I 1k;i1;m  ði  1Þði þ 1Þ I 1k;iþ1;m h io 2 þðjkÞ I 3i1;k;m  I 3iþ1;k;m þ ði þ 1ÞI 1k;i1;m  ði  1ÞI 1k;iþ1;m dljn ; h i bði; k; m; j; l; nÞ ¼ ði  1ÞI 2i1;k;m þ ði þ 1ÞI 2iþ1;k;m þ ði þ 1Þði  1ÞI 1i1;k;m  ði  1Þði þ 1ÞI 1iþ1;k;m djln ; aikm ¼ a1ikm  a2ikm ; bikm ¼ a1ikm  a2ikm ; (s (s ; k ¼ m  i; m þ i; ;  dikm ¼ 4 dikm ¼ 4 0; others; 0; Bmn wN ¼

b1;mn wNm2;n

þ b2;mn wmn N þ

k ¼ m  i; i  m; others;

b3;mn wNmþ2;n ;

b1;mn ¼ ða2m1 þ b2n Þð1 þ ðm þ 1Þ2 Þ  ða2mþ1 þ b2n Þð1 þ ðm  1Þ2 Þ; b2;mn ¼ ða2m1 þ b2n Þð1 þ ðm  3Þðm þ 1ÞÞ; b3;mn ¼ ða2mþ1 þ b2n Þð1 þ ðm þ 3Þðm  1ÞÞ;  K mn 2 ¼ kmn and

Amn 2;N ðwN ; wN Þ ¼

N X

wijN wkl N a2 ði; k; m; j; l; nÞ;

i;j;k;l¼1

a2 ði; k; m; j; l; nÞ ¼

16p4 ðl1 l2 Þ

7=2

ðl2 aði; k; m; j; l; nÞ  jlbði; k; m; j; l; nÞÞ;

1075

1076

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

for (8):

h i 2 2 2 2 aði; k; m; j; l; nÞ ¼ ði  1Þ I 1i1;m;k  ði þ 1Þ I 1iþ1;m;k  ði þ 1Þði  1Þ I 2i1;m;k þ ði  1Þði þ 1Þ I 2iþ1;m;k ajln ; h i bði; k; m; j; l; nÞ ¼ ði  1ÞI 4i1;k;m þ ði þ 1ÞI 4iþ1;k;m þ ði  1Þði þ 1ÞI 5i1;k;m  ði þ 1Þði  1ÞI 5iþ1;k;m bjln ;

for (9):

h i aði; k; m; j; l; nÞ ¼ ði  1Þ2 I 3i1;m;k  ði þ 1Þ2 I 3iþ1;m;k þ ði þ 1Þði  1Þ2 I 6i1;m;k  ði  1Þði þ 1Þ2 I 6iþ1;m;k dnjl ; h i bði; k; m; j; l; nÞ ¼ ði  1ÞI 7i1;k;m þ ði þ 1ÞI 7iþ1;k;m þ ði  1Þði þ 1ÞI 8i1;k;m  ði þ 1Þði  1ÞI 8iþ1;k;m bjln . I 4ikm ¼ ðk  1Þai;k1;m þ ðk þ 1Þai;kþ1;m þ ðk  1Þðk þ 1Þdk1;i;m  ðk þ 1Þðk  1Þdkþ1;i;m ; I 5ikm ¼ ðk  1Þdi;k1;m þ ðk þ 1Þdi;kþ1;m þ ðk þ 1Þðk  1Þbi;k1;m  ðk  1Þðk þ 1Þbi;kþ1;m ; I 6ikm ¼ bk1;m;i  bkþ1;m;i þ ðk þ 1Þdm;i;k1  ðk  1Þdm;i;kþ1 ;

I 7ikm ¼ ðk  1Þdm;i;k1 þ ðk þ 1Þdm;i;kþ1 þ ðk þ 1Þðk  1Þbk1;m;i  ðk  1Þðk þ 1Þbkþ1;m;i ; I 8ikm ¼ ðk  1Þbi;m;k1 þ ðk þ 1Þbi;m;kþ1 þ ðk þ 1Þðk  1Þdi;k1;m  ðk  1Þðk þ 1Þdi;kþ1;m ;   s cos pðm þ i þ kÞ  1 cos pðm  i  kÞ  1 1 þ aikm ¼ ; 4p mþiþk mik   s cos pðm þ i  kÞ  1 cos pðm  i þ kÞ  1 þ ; s ¼ l1 ; l2 . a2ikm ¼ 4p mþik miþk According to the obtained conclusions we will have mn mn mn K mn 1 wN ¼ A1;N ðwN ; wN Þ þ kB w ;

ð27Þ

kmn wmn N

ð28Þ

¼

Amn 2;N ðwN ; wN Þ;

m; n ¼ 1; 2; . . . ; N .

Lemma 2. For the coefficients of the expansion of functions W(x, y) and W(x, y), which satisfy (1), (7), (8) or (1), (7), (9) for fixed k, the following estimations hold:    2 2  wmn ¼ O m2 þ n2 ð1 þ m2 Þ1 ; wmn ¼ O m2 þ n2 . ð29Þ Proof. Since the series (16) converge, then wmn, wmn ! 0, when m, n ! 1. By virtue of the infinite system to (27) and (28), definitions a1(i, k, m, j, l, n), a2(i, k, m, j, l, n), K1, B, kmn , for fixed k we will get the formulae (29). h 3.3. On the convergence of the projective method Let us introduce the notation: mn emn  wmn 1 ¼ w N ;

emn 1

mn

emn 2

mn emn  wmn N ; 2 ¼ w mn

m; n ¼ 1; 2; . . . ; N ;

¼w ; ¼w ; m ¼ 1; 2; . . . ; 1; n ¼ N þ 1; N þ 2; . . . ; 1; m ¼ N þ 1; N þ 2; . . . ; 1; n ¼ 1; 2; . . . ; 1; m; n ¼ N þ 1; N þ 2; . . . ; 1;

ð30Þ

mn where femn 1 ; e2 g denotes the error of the solution of the truncated system relativity to the infinite one. Let us show, that for fixed k the eigenfunction of the truncated system (14) and (15) will converge to the eigenfunction of the infinite one (17) and (18).2

2

When we have Au  ku = 0, where A is a positive definite and self-adjointed operator it has been proved, that eigenelements may be obtained as limits of ‘‘approximate’’ eigenelements, obtained by the Galerkin method [10,14].

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

1077

mn mn mn Lemma 3. A sequence fwmn N ; wN ; m; n ¼ 1; 2; . . . ; 1ðwN ¼ wN ¼ 0; m ¼ 1; 2; . . . ; 1; n ¼ N þ 1; N þ 2; . . . ; mn 1; m ¼ N þ 1; N þ 2; . . . ; 1; n ¼ 1; 2; . . . ; 1; m; n ¼ N þ 1; N þ 2; . . . ; 1Þg, where fwmn N ; wN ; m; n ¼ 1; 2; . . . ; N g mn mn is an eigenfunction of (14) and (15) for fixed k converges to the sequence {w , w , m, n = 1, 2, . . . , 1}, satisfying (17) and (18). Moreover, for the simply supported plate 2

2 2 3 jemn i j 6 Cðm þ n Þ N

ð31Þ

and for the partially clamped plate 2 2 2 2 1 3 jemn i j 6 Cðm þ n Þ ð1 þ m Þ N

ð32Þ

(i = 1, 2; C = const > 0, which, obviously, depends on l1, l2). Proof. Obviously, (17) and (18) can be rewritten in the following form: mn mn mn K mn 1 w ¼ A1;N ðw; wÞ  kB w þ R1;N ðw; wÞ;

ð33Þ

kmn wmn

ð34Þ

¼

Amn 2;N ðw; wÞ



Rmn 2;N ðw; wÞ;

m; n ¼ 1; 2; . . . ;

where 1 X

Rmn s;N ðw; wÞ ¼

i;j;k;l¼N þ1

þ

1 X

wij wkl as ði; k; m; j; l; nÞ þ

N X

wij wkl as ði; k; m; j; l; nÞ þ   

i;j;k¼N þ1 l¼1

1 N X X

wij wkl as ði; k; m; j; l; nÞ þ    þ

i;j¼N þ1 k;l¼1

1 N X X

wij wkl as ði; k; m; j; l; nÞ þ    ;

s ¼ 1; 2

i¼N þ1 j;k;l¼1

and K 2 ¼ kmn in both cases, as we have just shown before. The number of combinations of sums in a definition of Rmn N is 15. Hence, in accordance with Lemma 1 and definitions as(i, k, m, j, l, n) we obtain 3 jRmn s;N ðw; wÞj 6 cN

ðc ¼ const > 0Þ.

From (14), (15) and (33), (34) by virtue of (30) we will have mn mn mn mn mn ðK mn 1 þ kB Þe1 ¼ A1;N ðe1 ; e2 Þ þ A1;N ðw; e2 Þ þ A1;N ðe1 ; wÞ þ R1;N ðw; wÞ; mn

mn

e e emn 2 ¼  A 2;N ðe1 ; e1 Þ  2 A 2;N ðw; e1 Þ  where

Rmn s;N ðw; wÞ

¼

e mn ðw; wÞ; R 2;N

m; n ¼ 1; 2; . . . ; 1;

Rmn s;N ;1 ðw; e2 Þ þ mn

ously given notations (30)),

mn Rmn s;N ;2 ðe1 ; e2 Þ þ Rs;N ;3 ðe1 ; wÞ e ¼ 1 Amn . Excluding e2 in A 2;N kmn 2;N

(Rmn s;N ;i

ð35Þ ð36Þ

are grouped sums according to the previ-

(35) and (36) we obtain

mn mn mn mn mn e mn e mn e mn e mn ðK mn 1 þ kB Þe1 ¼ A1;N ðe1 ; A 2;N ðe1 ; e1 ÞÞ  2A1;N ðe1 ; A 2;N ðw; e1 ÞÞ  A1;N ðe1 ; R 2;N ðw; wÞÞ  A1;N ðw; A 2;N ðe1 ; e1 ÞÞ mn

mn

mn mn mn e e  2Amn 1;N ðw; A 2;N ðw; e1 ÞÞ  A1;N ðw; R 2;N ðw; wÞÞ þ A1;N ðe1 ; wÞ þ R1;N ðw; wÞ.

Since we have shown that Rmn 1;N ! 0, N ! 1, then from the last equation it follows that e1 ! 0, when N ! 1. Therefore, the eigenfunction (wN, wN) tends to the eigenfunction (w, w) of the problem (17) and (18), which corresponds to given fixed value k. Using the previous conclusions we obtain the evaluations (31) and (32). h Now let us prove the theorem which allows to evaluate the approximate solution WN(x, y), WN(x, y). Assume dN,1(x, y) = W(x, y)  WN(x, y), dN,2 = W(x, y)  WN(x, y), kdk = max(x,y)2Djd(x, y)j. Theorem 1. The following inequalities hold:  m   m  m m o dN ;i    C1  6 C 1 ; o dN ;i  6 C 2 ; m ¼ 1; 2; 3; i ¼ 1; 2; kdN ;i kC 6 3 ;      3m 3m m m ox C N oy C N N  m   m  m m  o dN ;i      6 C 3 ;  o dN ;i  6 C 4 ; m ¼ 2; 3; C m ¼ const > 0; j ¼ 1; 2; 3; 4. j oxm1 oy  ox oy m1  N 3m N 3m C

C

ð37Þ

1078

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

Proof. First, we rewrite (16) in the form: W ðx; yÞ ¼

N X

wij xij ðx; yÞ þ r1N ðx; yÞ;

i;j¼1 N X

Wðx; yÞ ¼

ð38Þ ij

w uij ðx; yÞ þ

r2N ðx; yÞ;

i;j¼1

where 1 1 X X

r1N ðx; yÞ ¼

wij xij ðx; yÞ þ

i¼N þ1 j¼N þ1 1 1 X X

r2N ðx; yÞ ¼

1 N X X

wij xij ðx; yÞ þ

N 1 X X

i¼N þ1 j¼1

wij uij ðx; yÞ þ

i¼N þ1 j¼N þ1

1 N X X i¼N þ1 j¼1

wij xij ðx; yÞ;

ð39Þ

wij uij ðx; yÞ.

ð40Þ

i¼1 j¼N þ1

wij uij ðx; yÞ þ

N 1 X X i¼1 j¼N þ1

Consider one of the combinations of the sums in (39) and (40). By virtue of Lemmas 1 and 3 and definitions xij and uij the following inequality will be fulfilled:   X  c N 1 X   1 ij w xij ðx; yÞ 6 3 ðc1 ¼ const > 0Þ.   i¼1 j¼N þ1  N Analogous evaluations may be obtained, for uij and the other sums in (39) and (40). Therefore, from (39) and (40) we have c2 N3

kriN k 6

ði ¼ 1; 2; c2 ¼ const > 0Þ.

Then according to (38) dN ;1 ðx; yÞ ¼ W ðx; yÞ  W N ðx; yÞ ¼

N X

eij1 xij ðx; yÞ þ r1N ðx; yÞ;

i;j¼1

dN ;2 ðx; yÞ ¼ Wðx; yÞ  WN ðx; yÞ ¼

N X

eij2 uij ðx; yÞ þ r2N ðx; yÞ.

i;j¼1

By Lemma 3, and the obtained estimations for r1N and r2N we obtain kdN ;i kC 6

C 01 . N3

The necessary estimations (37) for partial derivatives follow from the definitions of the basis functions.

h

4. The iterative scheme Let k be such that k 2 [k0, k0 + Q] (k0 = k11 + e, e > 0). Divide [k0, k0 + Q] into M parts so that kr = kr1 + dr, r = 1, 2, . . . , M, kM = k0 + Q. In this section we propose an iterative scheme for the solution of the nonlinear system of algebraic equations with respect to wmn N , which are obtained from substitution of (15) into (14).For this purpose we can use the incremental loading, Runge–Kutta and Newton methods (these methods are developed, e.g., for nonlinear problems of mechanics in Odens monograph [17]). Let us assume dkr = h. Then we will have rþ1

r

wN ¼ wN þ

3 X i¼0

r

r

li zi ðwN ; kr Þ;

ð41Þ

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088 r

r

r

r

r

r

where z0 ¼ KðwN ; kr Þ, zi ¼ KðwN þ S i ; kr þ ai hÞ, S i ¼ r KðwN ; kr Þ r;cþ1 wN

¼

Pi1

r j¼0 bij zj ,

1079

r = 0, 1, . . . , M,

r r F 1 wN ðwN ; kr ÞF k ðwN ; kr Þ;

r;c

r;c

r;c

¼ wN F 1 c ¼ 0; 1; . . . ; cr  1; ð42Þ wN ðwN ; kr ÞF ðwN ; kr Þ; h i h i r;cr r ðwN ;kÞ mn ðwN ;kÞ 12 1N 21 22 2N , F k ðwN ; kÞ ¼ oF mnok , wN ¼ ðw11 where wN ¼ wN , F wN ðwN ; kÞ ¼ oFow m1 n1 N ; wN ; . . . ; wN ; wN ; wN ; . . . ; wN ; . . . ; N

T

wNN 1 ; wNN 2 ; . . . ; wNN N Þ , F(wN, k) = (F11, F12, . . . , F1N, F21, F22, . . . , F2N, . . . , FN1,FN2, . . . , FNN), Fij = Fij(wN, k), mn mn e F mn ¼ K mn 1 wN þ A1;N ðwN ; A 2;N ðwN ; wN ÞÞ þ kB wN . r;0

In order to take a good initial approximation wN for the solution of (42) for kr we use the Runge–Kutta formulae (see 17.5, [18]) owN ¼ F 1 wN ðwN ; kr ÞF k ðwN ; kr Þ. ok

ð43Þ

The inversion of the Jacobian F wN ðwN ; kÞ presents certain difficulties, which are connected with a requires of computer memory for the storage of the matrices, and the volume of computations for the inversion of Jacobian itself. If one would like to avoid these difficulties one can use hybrid iterative methods. For example, if we apply the nonlinear Seidel method and the incremental loading method together, instead of (43) we get owmn mn N ¼ F 1 ðwmn mn;wmn N ; kr ÞF mn;k ðwN ; kr Þ; N ok mn ¼ oF , then (41) will be written in the form: where F mn;wmn owmn N N

rþ1

3 X

r

mn wmn N ¼ wN þ

r

l i zi ;

r z0

r zi

¼K

rþ1 r rþ1 rþ1 r r m;n1 m;nþ1 12 w11 ; wmn ; . . . ; wNN N ; wN ; . . . ; wN N ; wN N ; kr

¼K

rþ1 w11 N

rþ1 rþ1 þ S i ; w12 N

þ

r S i ; . . . ; wmn N

rþ1

r þ S i ; wNm;nþ1 r

rþ1 r rþ1 rþ1 r r m;n1 m;nþ1 12 ; wmn ; . . . ; wNN w11 N ; wN ; . . . ; wN N ; wN N ; kr

K

¼

F 1 mn;wmn N  F mn;k

ð41 Þ

r ¼ 0; 1; . . . ; M; m; n ¼ 1; 2; . . . ; N ;

i¼0

! ;

r þ S i ; . . . ; wNN N r

r

!

þ S i ; kr þ ai h ;

!

rþ1 r rþ1 rþ1 r r m;n1 m;nþ1 12 mn NN w11 ; w ; . . . ; w ; w ; w ; . . . ; w N N N N N N rþ1 r rþ1 rþ1 r r m;n1 m;nþ1 12 w11 ; wmn ; . . . ; wNN N ; wN ; . . . ; wN N ; wN N

!

! .

Formula (42) with the application of the Seidel method for the external iterations and the Newton method for the internal iterations, is rewritten in the following way: ! r;cþ1 wmn N

r;cþ1 r;cþ1

r;c

r;cþ1

r;c

r;c

r;c

m;n1 1 11 12 mnþ1 ¼ wmn ; wmn ; . . . ; wNN N F mn;wmn wN ; wN ; . . . ; wN N ; wN N N

 F mn

r;cþ1 r;cþ1 r;cþ1 r;c r;c r;c m;n1 12 mnþ1 w11 ; wmn ; . . . ; wNN N ; wN ; . . . ; wN N ; wN N

! ;

c ¼ 0; 1; . . . ; cr  1; m; n ¼ 1; 2; . . . ; N .

ð42 Þ

r;cr

For the calculation of wmn N we use the formula (15) r;cr

e mn r;cr r;cr wmn N ¼  A N ðwN ; wN Þ;

m; n ¼ 1; 2; . . . ; N .

ð44Þ

1080

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

The iterative process is constructed so that first a solution for k = k0 = k11 + e(r = 0) is considered.For this 0;0 purpose an initial approximation wN is chosen.A certain number of iterations are needed to obtain a required 0;c0 0;c0 numerical accuracy0;cand to compute wN by (42) or (42*).By substituting the calculated values wN into (44) one 0 is able to calculate wN .Thus, we obtain an eigenfunction for k = k0.Furthermore, by (41) or (41*) we compute 0;c0 1;c1 1;0 1;0 1 0 1 wN , by taking into account of wN ¼ wN .Next wN ¼ wN , wN is an initial approximation for finding wN by for0;c0 1;c1 for computing wN mulae (42) or (42*) for k1 = k0 + h (it is also possible using wN as an initial approximation r;cr r;cr without applying (41) or (41*)). Continuing these steps we obtain an eigenfunction wN ; wN ; r ¼ 0; 1; . . . ; M. The convergence of the Newton method depends on the choice of the initial approximation. For appropriate choice of it the iterative process will converge. 0 It is known, that if F 1 w ðwÞ exists and is continuous, and if an initial approximation w is considered such 0 that, kF ðwÞk 6 g, then k

k

k w wk 6 gM

q2 1 k ; 1  q2

ð45Þ 2

M gL where kF 1 w ðwÞk 6 M, q ¼ 2 < 1, Mg 0

P1

k¼0 ðq 0

2k 1

Þ < r, Fw(w) satisfies the Lipschitz conditions with a con-

stant L in the ball U r ðwÞ ¼ fw 2 E ; kw  w k < rg, and EN is a set N-dimensional real vectors with a norm P P 1 N N 2 2 kwk ¼ j¼1 i¼1 wij . The Newton method has a quadratic convergence, if q < 1. (see, e.g., Samarsky N

and Gulin [18]). s;c r;cr r;cr r;cr r;cr r;cr r mn mn mn mn mn mn mn Supposing, that emn  wmn  wmn N ; m; n ¼ 1; 2; . . . ; N , e1 ¼ w ; e2 ¼ w ; m ¼ N þ 1; N þ 2; 1 ¼w N , e2 ¼ w s;cr

s;cr

mn . . . ; 1; n ¼ 1; 2; . . . ; 1; m ¼ 1; 2; . . . ; 1; n ¼ N þ 1; N þ 2; . . . ; 1; m; n ¼ N þ 1; N þ 2; . . . ; 1. emn are 1 ; e2 mn the errors of the projective–iterative method, where w denotes an exact solution for k = kr. Obviously, we have s;cr

r;cr

s;cr

s;cr

r;cr

s;cr

mn mn mn mn mn  wmn j emn 1 j ¼ jw N þ wN  wN j 6 je1 j þ j f1 j; mn mn mn mn mn j emn  wmn 2 j ¼ jw N þ wN  wN j 6 je2 j þ j f2 j.

s;cr

s;cr

mn mn mn Here femn 1 ; e2 g denotes the error of the projective method and ff1 ; f2 g is the error of the iterative method. Hence, from Lemmas 1–3, inequality (45) and the definition of the derivative F wN ðwN ; kÞ for appropriate initial r;0 2 approximation (wN ¼ Oððm2 þ n2 Þ Þ) we have the inequalities for the simply supported and partially clamped plate correspondingly:   cr r;cr C1 1 q2 1 mn þ j ei j 6 ; c 2 3 1  q2 r ðm2 þ n2 Þ N   cr r;cr C1 1 q2 1 j 6 þ ðC 1 ¼ const > 0Þ. j emn c i 2 3 1  q2 r ðm2 þ n2 Þ ð1 þ m2 Þ N r;cr

r;cr

r;cr

r;cr

Now let dN ;1 ðx; yÞ ¼ W ðx; yÞ  W N ðx; yÞ, dN ;2 ðx; yÞ ¼ Wðx; yÞ  WN ðx; yÞ. Theorem 2. The error functions of the projective–iterative method (41) and (42) ((41*) and (42*)) satisfy the estimations: s;cr

k dN ;i k 6 C 1



cr

1 q2 1 þ c N 3 1  q2 r

 ðC 1 ¼ const > 0Þ.

ð46Þ

Proof. By the previous expressions, notations, Theorem 1 and inequalities: r;cr 1 1 r;cr r;cr r;cr X X jwij  wijN jjxij ðx; yÞj; j dN ;2 ðx; yÞj 6 jwij  wijN jjuij ðx; yÞj j dN ;1 ðx; yÞj 6 i;j¼1

we obtain the necessary estimations (46).

i;j¼1

h

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

1081

5. Numerical examples The numerical examples demonstrate one of the possible approaches of the calculation of the initial approximations for the iterative scheme. This algorithm is based on Lemmas 1–3. Let us first consider the initial approximation, when k = k0 = k11 + e (e is a small positive number). By Lemmas 1–3 without restriction 0;0

0;0

ij of we can start the computations by choosing w11 N ¼ 1, wN ¼ 0; i; j 6¼ 1, (for instance, for simply supported plate the weight of the first coefficient is essential in the series (16), we can see it from Lemma 1 and (24) and (25). Further, note that for k = k12 + e the weight of the coefficient w12 becomes essential. So in that case 1;0

1;0

ij as the initial approximation one may take w12 N ¼ 1, wN ¼ 0; i 6¼ 1; j 6¼ 2. This way one may choose the initial approximations). Now let gi, i = 1, 2, 3, . . . be eigenvalues of (3) and (4) with g1 < g2 < g3 <   , where every eigenvalue can have some multiplicity (e.g., for the simply supported square plate, we suppose g1 = k11, g2 = k12 = k21, g3 = k22 end etc.) We will consider finding k within the segments [gi + e, gi+1], i = 1, 2, 3, . . . Now generalizing all aforesaid we can describe the following algorithm of choosing the initial approximations (for simplicity 0;0

ij suppose w11 N ¼ wN ):

I. for g1 < k 6 g2: ij 11 1. w11 N ¼ 1 (or wN ¼ 1), wN ¼ 0; i; j 6¼ 1; II. for g2 < k 6 g3: r0 ;cr

0

r0 ;cr

0

1. wijN ¼ wijN ( wijN is computed on the first step I, 1) ij 12 2. w12 N ¼ 1 (or wN ¼ 1), wN ¼ 0; i 6¼ 1; j 6¼ 2;

III. for g3 < k 6 g4: r1 ;cr

1

r1 ;cr

1

1. wijN ¼ wijN ( wijN is computed on the second step II, 1) 2.

wijN

3.

w22 N

¼

r0 ;cr 0 wijN

¼ 1, and etc.

r0 ;cr 0 ( wijN

wijN

is computed on the second step II, 2)

¼ 0; i; j 6¼ 2;

In the algorithm 1, 2, 3, . . . denote the branches of the solutions. Obviously, according to this approach we sort out solutions and follow the branches of the solutions. For the simply supported square plate a number of solutions on the segments [gi + e, gi+1] will be equal respectively to 2(i  1) + 3 (see Fig. 3). Let us demonstrate as an example the calculation the first tree branches of the solutions for the simply supported square plate l1 = l2 = 1. Let kW N k ¼ maxðxi ;y j Þ2GM jW N ðxi ; y j Þj, where GM = {(xi, yj), i, j = 1, 2, . . . , M} be a mesh on G. By the tables, the coefficients of the expansion WN, WN are presented. The maximum deflections are shown on Fig. 3. To decrease the volume of computations one can first apply the procedure for small N and then use the obtained data as initial approximations for the calculations with larger N. Supposing z ¼ r;cr 1

r;cr

max16i;j6N j wijN  wijN j. 0;0

Example I (scheme (42*)). k = k11 + e = 20.239 (e = 0.5); N = 9; the initial approximation: w11 N ¼ 0:288, 0;0

wijN ¼ 0; i; j 6¼ 1; c0 = 3 (a number of iterations); z = 2.7 · 106 (Tables 3 and 4). r;cr

r;cr

Here and below for simplicity: wij  wijN ; wij  wijN .

1082

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088 8

6

Maximum deflection

4

2 2

1

3

0

–2

–4

–6

–8 10

20

30

40

50 60 Parameter lambda

70

80

90

Fig. 3. The maximum of the deflection kWNk with N = 3, M = 10. Lines 1, 2, 3 are correspondingly the first three branches of the solutions.

Table 3 Coefficients wij (wij = 0 for even i, j) i, j

1

1 3 5 7 9

3

0.288392 0.000310 1.0 · 105 1.8 · 106 5.0 · 107

5

7 5

1.0 · 10 6.3 · 106 8.5 · 107 2.7 · 107 1.0 · 107

0.000310 5.6 · 105 6.3 · 106 1.5 · 106 4.9 · 107

9 6

1.8 · 10 1.5 · 106 2.7 · 107 1.1 · 107 4.7 · 107

5.0 · 107 4.9 · 107 1.0 · 107 4.7 · 107 2.7 · 106

7

9

Table 4 Coefficients wij (wij = 0 for even i, j) i, j

1

3

0.044788 0.001297 6.9 · 105 1.2 · 105 3.2 · 106

1 3 5 7 9

5 5

6.9 · 10 4.7 · 105 1.1 · 105 3.3 · 106 1.2 · 106

0.001297 0.000335 4.7 · 105 1.1 · 105 3.6 · 106

5

1.2 · 10 1.1 · 105 3.3 · 106 1.3 · 106 5.6 · 107

3.2 · 106 3.6 · 106 1.2 · 106 5.6 · 107 2.7 · 107

Example II (scheme (42*)). k = k12 + e = 49.848; N = 10. 1;0

1;0

1;0

1;0

1;0

ij 6 13 31 33 1. w11 (Tables 5 and 6), N ¼ 2:612, wN ¼ wN ¼ 0:153, wN ¼ 0:051, wN ¼ 0, i, j 5 1, 3; c1 = 9; z = 3.4 · 10 1;0

1;0

ij 6 (Tables 7 and 8). 2. w12 N ¼ 0:181 wN ¼ 0; i 6¼ 1; j 6¼ 2; c1 = 3; z = 1.2 · 10

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

1083

Table 5 Coefficients wij (wij = 0 for even i, j) i, j

1

3

5

7

9

1 3 5 7 9

2.635918 0.157102 0.010767 0.001404 0.000339

0.157102 0.051911 0.008015 0.001558 0.000444

0.010767 0.008015 0.002374 0.000637 0.000213

0.001404 0.001558 0.000637 0.000213 8.1 · 105

0.000339 0.000444 0.000213 8.1 · 105 3.4 · 105

Table 6 Coefficients wij (wij = 0 for even i, j) i, j

1

3

5

7

9

1 3 5 7 9

3.107522 0.016493 0.007382 0.001290 0.000323

0.016493 0.017486 0.006502 0.001769 0.000537

0.007382 0.006502 0.002712 0.000964 0.000359

0.001290 0.001769 0.000964 0.000425 0.000187

0.000323 0.000537 0.000359 0.000187 9.3 · 105

Table 7 Coefficients wij (wij = 0 for even i, j) i, j

2

4

6

8

10

1 3 5 7 9

0.179786 0.000302 3.3 · 106 8.1 · 107 2.7 · 107

0.000121 2.6 · 105 4.1 · 106 1.1 · 106 3.7 · 107

1.2 · 105 5.8 · 106 1.3 · 106 4.6 · 107 1.9 · 107

1.7 · 106 1.3 · 106 2.5 · 107 1.1 · 107 4.9 · 108

5.0 · 107 4.4 · 107 9.9 · 108 4.8 · 108 2.5 · 108

Table 8 Coefficients wij (wij = 0 for even i, j) i, j

1

3

0.041821 0.0023856 0.000143 2.5 · 105 6.9 · 106

1 3 5 7 9

5

0.002259 6.6 · 105 1.6 · 106 1.1 · 106 4.4 · 107

0.000301 0.000112 2.7 · 105 8.5 · 106 3.2 · 106

7

9 5

2.8 · 10 2.0 · 105 6.0 · 106 2.3 · 106 1.0 · 106

6.3 · 106 5.9 · 106 2.1 · 106 9.4 · 107 4.6 · 107

Table 9 Coefficients wij (wij = 0 for even i, j) i, j

1

3

5

7

9

1 3 5 7 9

4.213183 0.423933 0.044197 0.006110 0.001327

0.423932 0.173230 0.039539 0.008768 0.002322

0.044197 0.039539 0.016100 0.005116 0.001673

0.006110 0.008768 0.005116 0.002132 0.000832

0.001327 0.002322 0.001673 0.000832 0.000371

Example III (scheme (42*)). k = k22 + e = 79.457; N = 10; 2;0

2;0

2;0

2;0

2;0

ij 6 13 31 33 1. w11 (Tables 9 and 10), N ¼ 4:022, wN ¼ wN ¼ 0:378; wN ¼ 0:159, wN ¼ 0, i, j 5 1,3; c = 11; z = 8.6 · 10

1084

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

Table 10 Coefficients wij (wij = 0 for even i, j) i, j

1

3

5

7

9

1 3 5 7 9

6.926350 0.114573 0.012760 0.003671 0.001050

0.114572 0.030832 0.011324 0.005888 0.002139

0.012760 0.011324 0.007529 0.003954 0.001776

0.003671 0.005888 0.003954 0.002205 0.001131

0.001050 0.002139 0.001776 0.001131 0.000644

Table 11 Coefficients wij (wij = 0 for odd i and even j) i, j

2

4

6

8

10

1 3 5 7 9

1.550943 0.112668 0.009669 0.001147 0.000258

0.060875 0.019862 0.004531 0.001039 0.000309

0.006919 0.005297 0.001779 0.000567 0.000210

0.000998 0.001073 0.000470 0.000180 7.6 · 105

0.000261 0.000330 0.000166 7.2 · 105 3.4 · 105

Table 12 Coefficients wij (wij = 0 for even i, j) i, j

1

3

5

7

9

1 3 5 7 9

2.783091 0.068536 0.017119 0.003226 0.000821

0.165463 0.024702 0.000499 0.000381 8.9 · 105

0.014279 0.009098 0.002981 0.001036 0.000389

0.002108 0.002240 0.001069 0.000454 0.000199

0.000476 0.000678 0.000414 0.000209 0.000105

Table 13 Coefficients wij (wij = 0 for odd i, j) i, j

2

4

6

8

10

2 4 6 8 10

0.124582 4.3 · 105 1.9 · 105 1.3 · 106 3.9 · 107

4.3 · 105 1.5 · 105 4.4 · 106 8.3 · 107 3.0 · 107

1.9 · 105 4.4 · 106 2.2 · 106 5.9 · 107 2.4 · 107

1.3 · 106 8.3 · 107 5.9 · 107 1.2 · 107 6.0 · 108

3.9 · 107 3.0 · 107 2.4 · 107 6.0 · 108 5.1 · 108

Table 14 Coefficients wij (wij = 0 for even i, j) i, j

1

3

5

7

9

1 3 5 7 9

0.026812 0.003629 6.4 · 104 6.5 · 105 1.6 · 105

0.003629 0.0007115 6.9 · 105 2.3 · 106 1.2 · 107

6.4 · 104 6.9 · 105 7.2 · 105 1.8 · 105 6.4 · 106

6.5 · 105 2.3 · 106 1.8 · 105 5.1 · 106 2.1 · 106

1.6 · 105 1.2 · 107 6.4 · 106 2.1 · 106 9.6 · 107

2;0

2;0

2;0

2;0

ij 6 2. w12 (Tables 11 and 12), N ¼ 0:181; wN ¼ 0; i 6¼ 1; j 6¼ 2; c2 = 7; z = 9.8 · 10

3. w22 ¼ 0:131, wij ¼ 0; i; j 6¼ 24; c = 3; z = 2.7 · 106 (Tables 13 and 14).

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

1085

6. Neural network approximation of buckling loads In this last section a neural network procedure is proposed for the approximation of the buckling loads. The method uses examples, produced by means of the previously described exact mathematical algorithm, for the training of a back propagation neural network. The trained network is used in the sequel to approximate the nonlinear relation between plate quantities, like the plate dimensions l1, l2, and the corresponding buckling loads k = [k1, . . . , km]. Note that the same procedure can be used if these data are obtained from experiments ([22]). Here, a multi-layer back-propagation error driven neural network is used to learn the relation ðl1 ; l2 Þ ! k

ð47Þ

with all other data (material and plate thickness the same). (l1, l2) and the corresponding parameter vectors k from the previously calculated tables are used as training examples. In the production mode, after training, the nonlinear network reproduces the relation (l1, l2) ! k, i.e., for a given set of plate dimensions (different from the ones used in training) it gives a prediction for the buckling loads. The neural network model used the well-known Multilayer Perceptron (MLP). It is the most widely used neural network model for function approximation with numerous successful applications in almost every scientific and engineering domain. The most attractive feature of the MLP is that it exhibits excellent interpolation capabilities (even when trained with sparse datasets) which make it an ideal solution for data-driven inverse modelling problems. The MLP model (also known as backpropagation neural network) is a feedforward neural network with one or more hidden layer containing units (called hidden units) with nonlinear activation function (usually of sigmoid type). In our experiments, to implement a mapping from a d-dimensional input space

Fig. 4. Training of the neural network for the buckling problem.

1086

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088 70

60

50

40

30

20

10

0

0

10

20

30

40

50

60

70

Fig. 5. Accuracy of neural network prediction for the buckling loads.

Table 15 e k 21 is prediction, Er (%) denotes the relativity error in percent e l1 k 21 (Nnet) l2

k21

Er (%)

1 1 1 1 1 2.5 2.5 2.5 2.5 3.5 3.5 3.5 3.5 4.5 4.5 5 5 5

66.5526 49.711 42.8261 41.1391 40.2798 44.634 10.6484 8.2355 7.40405 42.2283 8.61216 4.41792 4.15707 41.1924 4.53335 41.0008 2.99883 2.6621

0.120206 5.27127 6.2281 3.02387 1.83981 3.17779 1.22818 4.13697 5.75833 3.1844 4.94831 7.58468 1.444 1.56342 0.133403 1.50313 5.43448 8.00853

1 1.5 2.5 3.5 5 1 2.5 3.5 4.5 1 2.5 4.5 5 1 3.5 1 4.5 5

66.6326 52.3314 45.4934 42.3831 41.0209 43.2156 10.7792 8.5762 6.9777 40.8836 8.1860 4.7530 4.2171 40.5484 4.5273 40.3845 3.1618 2.8753

to a m-dimensional output space, we have used an MLP with d inputs, m outputs and one hidden layer with H hidden units with the hyperbolic tangent sigmoid activation function. The MLP model can be trained to implement the desired mapping from a d-dimensional to an m-dimensional domain by using a training set that contains N examples of the mapping, i.e., pairs of the form (l, k). For effectiveness we use different neural network for each buckling load ki, i.e., we consider d = 2 inputs (l1 and l2) and m = 1 outputs. Once a training set is available, the training process is actually an optimization

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

1087

procedure that adjusts the network parameters (weights and biases) to minimize a quadratic error function in the least squares sense. In our method the Levenberg–Marquadt method has been used for the minimization of the error function available from the Matlab Neural Network toolbox. This training method has been found to be the most effective among several tested training techniques and achieved to provide near zero minima of the error function even in the case of networks with small number H of hidden units. Finally, for every problem examined, after the completion of training, the prediction accuracy of the constructed network was assessed by using a separate test set of cases (different from the training set). Let us present a concrete numerical example: a 2–8–1 neural network is trained in the training data. Training up to an accuracy equal to 0.0001 requires 808 epochs (cycles of approximations, where each cycle involves parameter updating based on the whole set of available training data). The convergence of training is shown in Fig. 4. The accuracy of predictions for the learning data, i.e., k21 exact us. Neural network prediction e k 21 is shown in Fig. 5 and Table 15. Finally we have tested the network on different combinations of l1, l2. Some of the values k with the corresponding error are shown in Table 15. Acknowledgement The support of the Greek Ministry of Education through a research grant of the National Fellowships Foundation (IKY) is gratefully acknowledged. References [1] Allgower EL, Georg K. Numerical continuation. An introduction. Berlin: Springer; 1990. [2] Bauer L, Reiss E. Nonlinear buckling of rectangular plates. SIAM 1965;13:603–27. [3] Berger MS. On von Karmans equations and the buckling of a thin elastic plate, I, The clamped plate. Commun Pure Appl Math 1967;20:687–719. [4] Chien CS, Chang SL, Mei Z. Tracing the buckling of a rectangular plate with the block GMRES method. J Comput Appl Math 2001;136:199–218. [5] Chien CS, Kuo YJ, Mei Z. Bifurcations of the von Karman equations with Robin boundary conditions. Comput Math Appl 1999;38:85–112. [6] Ciarlet PG, Destuynger P. A justification of a nonlinear model in plate theory. Comput Methods Appl Mech Eng 1979;9:227–58. [7] Ciarlet P, Rabier P. Les equations de von Karman. Berlin, Heidelberg, New York: Springer-Verlag; 1980. [8] Caloz G, Rappaz J. Numerical analysis for nonlinear and bifurcation problems. Handbook of numerical analysis, vol. V. Amsterdam: North Holland; 1997. p. 487–637. [9] Crouzeix M, Rappaz J. On numerical approximation in bifurcation theory. Recherches en Mathe´matiques Applique¨es, 13. Masson, Paris. Berlin: Springer-Verlag; 1990. [10] Demidovich BP, Maron IA. The basis of computational mathematics, Gosdr. Izd. Physic.-Mathem. Liter., Moscow, 1960 (in Russian). [11] Dym CL, Shames IH. Solid mechanics, a variational approach. McGraw-Hill; 1973. p. 280–349. [12] Holder EJ, Schaeffer DG. Boundary conditions and mode jumping in the von Karmans equations. SIAM J Math Anal 1984;15:446–58. [13] Iyengar NGR. Structural stability of columns and plates. Ellis Horwood Limited; 1988. [14] Kantorovich LV, Krilov VI. Approximate methods of high analysis. Izd. Physic.-Mathem. Liter., Moscow, 1967 (in Russian). [15] Lions JL. Quelques Me´thodes de Re´solution des Proble`mes aux Limites non Line´aires. Paris, Dunod Gauthier-villars, 1969. [16] Muradova A. A design of some non-linear eigenvalue problems. Bulletin of TICMI (Tbilisi International Center of Mathematics and Informatics), 2002;6:24–7. [17] Oden JT. Finite elements of nonlinear continua. New York: McGraw-Hill Book Company; 1972. [18] Samarsky AA, Gulin AV. Numerical methods. Moskow: Nauka; 1989 (in Russian). [20] Schaeffer DG, Golubitsky M. Boundary conditions and mode jumping in the buckling of a rectangular plate. Commun Math Phys 1979;69:209–36. [21] Vashakmadze TS. The theory of anizotropic elastic plates. Dordrecht, Boston, London: Kluwer Academic Publishes; 1999. [22] Waszczyszyn Z, Ziemian´ski L. Neural networks in mechanics of structures and materials—new results and prospects of applications. Comput Struct 2001;79:2261–76. [23] Weinberg MM, Trenogin VA. The bifurcation theory of nonlinear equations. Moskow: Nauka; 1969 (in Russian). [24] Stavroulakis GE, Abdalla KM. A systematic neural network classificator in mechanics. An application in semi-rigid steel joints. Int J Eng Anal Des 1994;1:279–92.

1088

A.D. Muradova, G.E. Stavroulakis / Communications in Nonlinear Science and Numerical Simulation 12 (2007) 1068–1088

[25] Stavroulakis GE, Antes H. Nondestructive static unilateral crack identification. A BEM—neural network approach. Comput Mech 1997;20(5):439–51. [26] Stavroulakis GE. Inverse and identification problems in mechanics. Dordrecht, Boston, London: Kluwer Academic Publishers; 2000. [27] Waszczyszyn Z, editorNeural networks in the analysis and design of structures. CISM courses and lectures No. 404. Wien-New York: Springer; 1999. p. 1–52. and p. 161–96. (chapter 1, Fundamentals of artificial neural network, chapter 4, The neural network approach in plasticity and fracture mechanics, together with P.D. Panagiotopoulos). [28] Rojas R. Neural networks—a systematic introduction. Berlin, Heidelberg: Springer; 1996. [29] Haykin S. Neural networks—a comprehensive foundation. Upper Saddle River, NJ: Prentice-Hall; 1999.