GWW4r*aZ3. (12.83)

=1 + 4 Р = 2(М. + Л) Y = +

J = 1хЦу - уЦху V4 = Ixx + lyy и Т. Д.

Если различные члены, например а , в (12.79) определены в каждой точке сетки, можно легко осуществить дискретизацию (12.79) с центральными разностями.

Как отмечено в § 12.2, проще начать с численного расчета параметров xi и т. д., определенных выражениями (12.30) и (12.31). Различные члены в (12.79) могут быть выражены через х и т. д. Члены а , р , y/J могут быть получены из (12.12):

GTT- = .

(ШТ1 = Щ. (12.80)

Gww = = ibi.

где якобиан обратного преобразования/ = ху - Обозначения GTT и т. д. введены в соответствии с программой LAGEN (рис. 12.9).

Все члены правой части уравнений (12.80) могут быть получены из (12.30). Исходя из двумерного представления уравнений (12.7) для производной от 1ху можно получить

хУ1г\ + хУцц Ul + lxxri (12.81)

Получив аналогичные выражения для 1уу и сделав некоторые преобразования, члены (Vl/J) и (Vt) ) можно представить в виде

J J J

+ .2ZBfnZif!mL. (12.82)

DELETE -(ir I 1

3 С LAGEN APPLIES THE FINITE DIFFERENCE METHOD TO LAPLACES

4 С EQUATION IN GENERALISED COORDINATES ON A MODIFIED POLAR 6Ш

5 С THE DISCRETISED EQUATION IS SOLVED BY SOR

7 DIMENSION XG(23,23),YG(23,23),GWW(21,21),GWT(21,21),CTT(21,2l)r

8 1DELZI(21.21),DELET(21,21).PHI(21,21),PHIX(21,21)

9 COMMON /GRIDP/XG,YG,PHIX,PHI

10 COMMON /TRAPP/GWW.GWT.GTT.DELZI.DELtT

11 С

12 OPEN(1,FILE=LAGEN.DAT)

13 OPEN(6,FILE LAGEN.OUT)

14 READ(1,1)JMAX.KMAX,NMAX,lEX.IPR

15 READ(1,2)RW,RX,RY,RZ,THEB,THEN,EPS,0H

16 READ(1,2)BXI,BXXI,CXI,CXXI

17 1 FORMAT(815)

18 2 FORMAT(8E10.3)

19 С

20 WRITE(6,3)

21 WRITE(6,4)JMAX,KMAX.NMAX,lEX,EPS,OK

22 VRITE(6,5)RV,RX,RY,RZ,THEB,THEN

23 3 FORMAT С LAPLACE EQUATION BY GEM. COORD. FDH*, )

24 4 FORMATC JMAX*,12, KMAX ,I2, ИМАХ М5/ IEX M2,

25 1 EPS ,E10.3, OM=,F5.3)

26 5 FORMATC RW ,F5.3, RX ,F5.3/ RY-M5.3, RZ *,F5.3r

27 15X, THEB=,F5.1, THEM*,F5.1, )

28 ARM KMAX - 2

29 AJM = JMAX - 2

30 С

31 С SET UP GRID

32 С

33 CALL GRID(JMAX,KMAX,THED,THEN,RW,RX,RY,RZ)

34 С

35 С SET BOUNDARY VALUES OF PHI

36 с

37 DO 8 J 1,JMAX

38 PHKJ.l) 0.

39 РНК J, KMAX) PHIX(J,KMAX)

40 8 CONTINUE

41 DO 9 К = 1,KMAX

42 PHI(1,K) PHIX(1,K)

43 PHI(JMAX,K) PHIX(JMAX,K)

44 9 CONTINUE

45 С

46 С SET GRID RELATED PARAMETERS

47 С

48 CALL TRAPA(JMAX,KMAX)

49 С

50 С ITERATE USING SOR

51 С

52 CALL ITER(JMAX,KMAX,NMAX,N,OM,EPS)

53 С

54 С COMPARE SOLUTION WITH EXACT

55 С

56 10 SUM = 0.

57 DO 15 К 1,KMAX

58 VRITE(6,11)K

59 11 FORMAT(/, K M2)

60 DO 12 J = 1,JMAX

61 DIF = PHI(J,K) - PHIX(J,K)

62 SUM = SUM + DIF*DIF

63 12 CONTINUE

64 VRITE(6,13)(PHI(J.K),J=1,JMAX)

65 WRITE(6,14)(PHIX(J,K),J=1,JHAX)

66 13 FORMATC PHI=,10F7.4)

67 14 FORMATC PHX=,10F7.4)

68 15 CONTINUE

69 RMS = SQRT(SUM/AJM/AKN)

70 VRITE(6,16)N,RMS

71 16 FORMAT(/, CONVERGED AFTER ,13, STEPS, RMS \R12.5)

72 17 CONTINUE

73 STOP

74 END

Рис. 12.9. Распечатка программы LAGEN.

2 subroutine GRID(jmax.киах,тнев,THEN,rw,RX,ry,rz)

4 с set the augmented grid, initial and exact phi

6 dimension xg(23.23),YG(23.23),phix(21,21),phi(21,21)

7 common /gridp/xg,YG,phix,phi

8 jmap = jmax - 1

9 kmap = kmax - 1

10 ajm = jmap

11 akm kmap

12 drwx = (rx - rv)/ajm

13 drzy = (ry - rz>/akm

14 dth = (then-theb)/akm

15 pi = 3.1415927

16 kpp = kmax + 2

17 гр = jmax + 2

18 с

19 с SE; XG, YG. exact and initial phi

20 с

21 do 7 к l.kpp

22 ak = к - 2

23 thk = (theb + ak*DT:n*ri/180.

24 ck = cos(thk)

25 sk = sin(thk)

26 dr = drwx + (drzy - drwx)*АК/ЛКМ

27 rwz = rw + (rz - RW)*ak/akm

28 do 6 j = 1,jpp

29 aj = j - 2

30 r = rwz + aj*dr

31 xg(j,k) = r*ck

32 yg(j,k) = r*sk

33 if(k .eq. 1 .or. к .eq. kpp)goto g

34 if(j .eq. 1 .or. j .eq. jpp)goto 6

35 jm = j-1

36 km к - 1

37 phix(jm.km) = sk/r

38 phkjm.km) = phix(jm,km)

39 6 continue

40 7 continue

41 return

42 end

Рис. 12.10. Распечатка подпрограммы GRID.

Все члены, стоящие в правой части (12.82), (12.83), можно выразить через (12.30), (12.31) и (12.80). После определения QTT и т. д. через (12.80) - (12.83) в каждом узле сетки уравнение (12.79) может быть при помощи трехточечных центральных разностей приведено к виду

0.5 [(DELZI.),+i (DELZI.</>)/ ,] -

- 0.5 [(DELET.<)/, - (DELET.</>)/, +

+ [(QTT.</>)/i - 2 (QTT.<)/ + (QTT <)/+ ,] +

f 0.25[(QWT.)y+ (GWT )/i, + (QWT <)/ , -