home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / prgramer / adaptor / examples / pde / expl21.f next >
Text File  |  1993-06-28  |  24KB  |  714 lines
  1.       PROGRAM EXPL2
  2. C***********************************************************************
  3. C                                                                      *
  4. C AUTHOR: YVONNE LUH, GMD, I1.HR                                       *
  5. C                                                                      *
  6. C         POSTFACH 1319                                                *
  7. C         W-5205 ST. AUGUSTIN                                          *
  8. C                                                                      *
  9. C***********************************************************************
  10. C***********************************************************************
  11. C     fortran 90 - Version
  12. C     28/4/92
  13. C.......................................................................
  14. C     EXPLICIT ALGORITHM: STANDARD 2ND ORDER DISCRETIZATION
  15. C     WAVE EQUATION (2 SPACE DIMENSIONS)
  16. C
  17. C     INITIAL BOUNDARY VALUE PROBLEM
  18. C     DDU/DDT=DDU/DDX+DDU/DDY
  19. C     U(X,Y,0)     = F1(X,Y)
  20. C     DU/DT(X,Y,O) = F2(X,Y)
  21. C     U(X,Y,T)     = G(X,Y,T)    (X,Y) AT BOUNDARY OF (0,1)X(0,1)
  22. C.......................................................................
  23. C     YOU CAN TRY DIFFERENT K,H
  24. C     K    STEP IN T-DIRECTION
  25. C     H    STEP IN X- AND Y-DIRECTION
  26. C.......................................................................
  27. C     DIFFERENT EXAMPLES CAN BE CALCULATED
  28. C               COMMON /NEXAMP /NEXAMP
  29. C                1      DSIN(PI*X)*DSIN(PI*Y)
  30. C                            * DCOS(DSQRT(2.0)*PI*T)
  31. C                            * 1000.0D0
  32. C                2      DSIN(20.0*PI*X)*DSIN(20.0*PI*Y)
  33. C                            * DCOS(20.0D0*DSQRT(2.0D0)*PI*T)
  34. C                3      DEXP(X+Y)*DEXP(DSQRT(2.0D0)*T)
  35. C                            * 100000.0D0
  36. C                4      DEXP( 3.0D0*(X+Y))
  37. C                            * DEXP( 3.0D0*DSQRT(2.0D0)*T)
  38. C                5      X**2 + Y**2 + 2.0D0*T**2
  39. C                6      0.0
  40. C     DIFFERENT KINDS OF ERRORS CAN BE CALCULATED
  41. C               COMMON /NRELER  /NRELER
  42. C                1      RELATIVE ERROR    =MAX((SOL - U)/SOL)
  43. C                2      ABSOLUTE ERROR    =MAX (SOL - U)
  44. C***********************************************************************
  45. C***********************************************************************
  46. c
  47.       double precision   u0(:,:), u1(:,:), u2(:,:), tx(:,:), ty(:,:)
  48.       double precision   hu1 (:,:)
  49.       double precision   hu2 (:,:)
  50.       DOUBLE PRECISION   TEND,T,K,H,P
  51.       DOUBLE PRECISION   mxresul, l2resul
  52.       INTEGER            NP,N
  53.       INTEGER            NEXAMP,NRELER
  54.       INTEGER            I, J
  55. c
  56. C.......................................................................
  57. c
  58. C     COMMON BLOCKS
  59. C     NEXAMP             NUMBER OF EXAMPLE
  60. C     NRELER             RELATIVE OR ABSOLUTE ERROR IS CALCULATED
  61. c
  62.       COMMON /NEXAMP /NEXAMP
  63.       COMMON /NRELER /NRELER
  64. c
  65. C***********************************************************************
  66. C***********************************************************************
  67. C     INPUT OF PROBLEM PARAMETERS
  68. C***********************************************************************
  69. C***********************************************************************
  70. C
  71.       READ(5,*) K
  72.       READ(5,*) NP
  73.       READ(5,*) TEND
  74.       READ(5,*) NEXAMP
  75.       READ(5,*) NRELER
  76. C
  77. C***********************************************************************
  78. C***********************************************************************
  79. C     declaration of allocatable arrays
  80. C***********************************************************************
  81. C***********************************************************************
  82. C
  83.       allocate  (u0(0:np,0:np), u1(0:np,0:np), u2(0:np,0:np),
  84.      +           tx(0:np,0:np), ty(0:np,0:np))
  85.       allocate (hu1(0:np,0:np))
  86.       allocate (hu2(0:np,0:np))
  87. C
  88. C***********************************************************************
  89. C***********************************************************************
  90. C     COMPUTATION OF PARAMETERS
  91. C***********************************************************************
  92. C***********************************************************************
  93. C
  94.       N=NP-1
  95.       H=1.0D0/DBLE(NP)
  96.       P=K/H
  97. C
  98. C***********************************************************************
  99. C***********************************************************************
  100. C     calculation of arrays tx and ty
  101. C***********************************************************************
  102. C***********************************************************************
  103. C
  104. !HPF$ INDEPENDENT, LOCAL_ACCESS
  105.       do j = 0, np 
  106. !HPF$ INDEPENDENT, LOCAL_ACCESS
  107.         do i = 0, np
  108.           tx(i,j) = dble (i)
  109.           ty(i,j) = dble (j)
  110.         end do
  111.       end do
  112.         tx = tx * h
  113.         ty = ty * h
  114. c
  115. C***********************************************************************
  116. C***********************************************************************
  117. C     OUTPUT OF PARAMETERS
  118. C***********************************************************************
  119. C***********************************************************************
  120. C
  121.       WRITE(6,*)'*****************************************************'
  122.       WRITE(6,*)'EXPLICIT DISCRETIZATION OF WAVE EQUATION (2D)'
  123.       WRITE(6,*)'*****************************************************'
  124.       WRITE(6,*) ' COMPUTED EXAMPLE: '
  125.       IF      (NEXAMP.EQ.1) THEN
  126.            WRITE(6,*) ' SIN(PI*X)*SIN(PI*Y)*COS(1.414*PI*T) *1000'
  127.       ELSE IF (NEXAMP.EQ.2) THEN
  128.            WRITE(6,*) ' SIN(20*PI*X)*SIN(20*PI*Y)*COS(20*1.414*PI*T)'
  129.       ELSE IF (NEXAMP.EQ.3) THEN
  130.            WRITE(6,*) ' EXP(X+Y)*EXP(1.414*T) * 100000'
  131.       ELSE IF (NEXAMP.EQ.4) THEN
  132.            WRITE(6,*) ' EXP( 3*(X+Y))*EXP( 3*1.414*T)'
  133.       ELSE IF (NEXAMP.EQ.5) THEN
  134.            WRITE(6,*) ' X**2 + Y**2 +2*T**2'
  135.       ELSE IF (NEXAMP.EQ.6) THEN
  136.            WRITE(6,*) ' 0.0 '
  137.       ENDIF
  138.       IF      (Nreler.EQ.1) THEN
  139.            WRITE(6,*) ' relative error is calculated' 
  140.       ELSE IF (Nreler.EQ.2) THEN
  141.            WRITE(6,*) ' absolute error is calculated' 
  142.       ENDIF
  143.       WRITE(6,*) '****************************************************'
  144.       WRITE(6,9000) ' TIME -STEP K = ', K
  145.       WRITE(6,9000) ' SPACE-STEP H = ', H
  146.       WRITE(6,9000) ' RATIO P=K/H  = ', P
  147. C
  148.       IF (P*P.GT.0.5D0) THEN
  149.          WRITE(6,*) ' '
  150.          WRITE(6,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
  151.          WRITE(6,*) ' CFL CONDITION (P*P =< 0.5) NOT FULFILLED'
  152.          WRITE(6,*) ' THIS WILL CAUSE LARGE ERRORS IN COMPUTED SOLUTION'
  153.          WRITE(6,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
  154.          WRITE(6,*) ' '
  155.       ENDIF
  156. C
  157.       WRITE(6,9001) ' NUMBER OF GRID POINTS PER LINE = ',NP
  158.       WRITE(6,9000) ' MAX. T       = ', TEND
  159.  9000 FORMAT(1X,A17,F20.10)
  160.  9001 FORMAT(1X,A35,I5)
  161.       WRITE(6,*) '*****************************************************'
  162. C
  163. C***********************************************************************
  164. C***********************************************************************
  165. C     COMPUTATION OF INITIAL VALUES IN FIRST TWO TIME STEPS
  166. C***********************************************************************
  167. C***********************************************************************
  168. C
  169.       call f1(u0,tx,ty,np)
  170.       call sol(u1,tx,ty,np,k)
  171.       T=K
  172. C
  173. C***********************************************************************
  174. C***********************************************************************
  175. C     TEST OF INITIAL VALUES AT TIME T=K
  176. C***********************************************************************
  177. C***********************************************************************
  178. C
  179.       call difmx(u1,tx,ty,np,k,mxresul)
  180.       call difl2(u1,tx,ty,np,k,l2resul)
  181.       WRITE(6,9020) ' T = ', ' difmx ', ' difl2 '
  182.       WRITE(6,9021) t, mxresul, l2resul
  183.  9020 FORMAT (1X,A8,2x,a10,2x,a10)
  184.  9021 FORMAT (1X,F8.4,2x,D10.4,2x,D10.4)
  185. C
  186.       T=K+K
  187. C
  188. C***********************************************************************
  189. C***********************************************************************
  190. C
  191.    40 CONTINUE
  192. C
  193. C***********************************************************************
  194. C***********************************************************************
  195. C     COMPUTATION OF SOLUTION
  196. C***********************************************************************
  197. C***********************************************************************
  198. C
  199. C     IN INTERIOR
  200. C
  201. C     overlap U1(:,1:1) with HU1
  202. C
  203. C     n = np - 1
  204.  
  205.       hu1(0:np,1:n) = u1(0:np,2:n+1)
  206.       hu2(0:np,1:n) = u1(0:np,0:n-1)
  207.  
  208. !HPF$ INDEPENDENT, LOCAL_ACCESS
  209.       do j=1,n
  210. !HPF$ INDEPENDENT, LOCAL_ACCESS
  211.         do i=1,n
  212.          U2(i,j) = P*P* ( U1(i+1,j) + U1(i-1,j) 
  213.      +        + Hu1(i,j) + HU2(i,j)
  214.      +        - 4.0D0 * U1(i,j) )
  215.      +        + 2.0D0 * U1(i,j) - U0(i,j)
  216.        end do
  217.       end do
  218. C
  219. C.......................................................................
  220. C
  221. C     ON BOUNDARY
  222. C
  223.       call g(u2,tx,ty,np,t)
  224. C
  225. C***********************************************************************
  226. C***********************************************************************
  227. C     COMPUTATION AND OUTPUT OF ERROR OF THE COMPUTED SOLUTION
  228. C***********************************************************************
  229. C***********************************************************************
  230. C
  231.       call difmx(u2,tx,ty,np,t,mxresul)
  232.       call difl2(u2,tx,ty,np,t,l2resul)
  233.       WRITE(6,9021) t, mxresul, l2resul
  234. C
  235. C***********************************************************************
  236. C***********************************************************************
  237. C     STORAGE OF NEW SOLUTION
  238. C***********************************************************************
  239. C***********************************************************************
  240. C
  241.       U0 = U1
  242.       U1 = U2
  243. C
  244. C***********************************************************************
  245. C***********************************************************************
  246. C     NEXT TIME STEP
  247. C***********************************************************************
  248. C***********************************************************************
  249. C
  250.       T=T+K
  251.       IF (T.LT.TEND) GOTO 40
  252. C
  253. C***********************************************************************
  254. C***********************************************************************
  255. C     deallocation
  256. C***********************************************************************
  257. C***********************************************************************
  258. c
  259.       deallocate (hu2, hu1)
  260.       deallocate (ty, tx, u2, u1, u0)
  261. c 
  262. C***********************************************************************
  263. C***********************************************************************
  264. C
  265.       END
  266. C
  267. C****************************************************************
  268. C****************************************************************
  269. C****************************************************************
  270. C****************************************************************
  271. C****************************************************************
  272. C
  273. C
  274.       subroutine DIFMX(U,tx,ty,Np,T,mxresul)
  275. C
  276. C     COMPUTES THE MAXIMUM NORM OF THE DIFFERENCE BETWEEN
  277. C     SOL (=SOLUTION OF THE BVP) AND THE VALUES IN U
  278. C
  279. C     DEPENDING ON NRELER  THE RELATIVE ERROR (NRELER=1)
  280. C     OR THE ABSOLUTE ERROR (NRELER=2) IS COMPUTED
  281. C
  282.       integer          np, n
  283.       DOUBLE PRECISION U(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
  284.       DOUBLE PRECISION mxresul
  285.       DOUBLE PRECISION uhelp  (0:np,0:np)
  286.       DOUBLE PRECISION uhelp1 (0:np,0:np)
  287.       INTEGER          t
  288.       INTEGER          NRELER
  289.       COMMON /NRELER /NRELER
  290. C
  291.       N = Np-1
  292. C
  293.       call sol(uhelp,tx,ty,np,t)
  294. c      
  295.       uhelp1 = uhelp
  296.  
  297.       if (nreler.eq.1) then
  298.           where ( uhelp(1:n,1:n) .gt. 1.0d-10 )      
  299.               uhelp1(1:n,1:n) = abs( (uhelp(1:n,1:n)-u(1:n,1:n))
  300.      +                              / uhelp(1:n,1:n) )
  301.           elsewhere
  302.                uhelp1(1:n,1:n) = abs( uhelp(1:n,1:n) - u(1:n,1:n) )
  303.           endwhere
  304.       else if (nreler.eq.2) then
  305.                uhelp1(1:n,1:n) = abs( uhelp(1:n,1:n) - u(1:n,1:n) )
  306.       endif
  307. c
  308.       mxresul = maxval (uhelp1(1:n,1:n))
  309.     
  310.       END
  311. C
  312. C
  313. C****************************************************************
  314. C****************************************************************
  315. C****************************************************************
  316. C****************************************************************
  317. C
  318. C
  319.       subroutine DIFl2(U,tx,ty,Np,T,l2resul)
  320. C
  321. C     COMPUTES THE L2 NORM OF THE DIFFERENCE BETWEEN
  322. C     SOL (=SOLUTION OF THE BVP) AND THE VALUES IN U
  323. C
  324. C     DEPENDING ON NRELER  THE RELATIVE ERROR (NRELER=1)
  325. C     OR THE ABSOLUTE ERROR (NRELER=2) IS COMPUTED
  326. C
  327.       integer          np, n
  328.       DOUBLE PRECISION U(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
  329.       DOUBLE PRECISION l2resul
  330.       DOUBLE PRECISION uhelp(0:np,0:np)
  331.       DOUBLE PRECISION uhelp1(0:np,0:np)
  332.       INTEGER          t
  333.       INTEGER          NRELER
  334.       COMMON /NRELER /NRELER
  335. C
  336.       N = Np-1
  337. C
  338.       call sol(uhelp,tx,ty,np,t)
  339. c
  340.       uhelp1 = uhelp
  341.       if (nreler.eq.1) then
  342.           where ( uhelp(1:n,1:n) .gt. 1.0d-10 )
  343.               uhelp1(1:n,1:n) = abs( (uhelp(1:n,1:n)-u(1:n,1:n))
  344.      +                              / uhelp(1:n,1:n) )
  345.           elsewhere
  346.                uhelp1(1:n,1:n) = abs( uhelp(1:n,1:n) - u(1:n,1:n) )
  347.           endwhere
  348.       else if (nreler.eq.2) then
  349.                uhelp1(1:n,1:n) = abs( uhelp(1:n,1:n) - u(1:n,1:n) )
  350.       endif
  351. c
  352.       uhelp1(1:n,1:n) = uhelp1(1:n,1:n) ** 2
  353.       l2resul = sum ( uhelp1(1:n,1:n) )
  354. c
  355.       l2resul = dsqrt(l2resul) / dble(n)
  356. c
  357.       END
  358. C
  359. C
  360. C****************************************************************
  361. C****************************************************************
  362. C****************************************************************
  363. C****************************************************************
  364. C****************************************************************
  365. C
  366. C
  367.       subroutine SOL(u,tx,ty,np,T)
  368. C
  369. C**********************************************************************
  370. C
  371. C     CONTINUOUS SOLUTION OF THE DIFFERENTIAL PROBLEM
  372. C
  373. C**********************************************************************
  374. C
  375.       integer          np
  376.       DOUBLE PRECISION u(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
  377.       DOUBLE PRECISION t
  378.       DOUBLE PRECISION PI
  379.       INTEGER          NEXAMP
  380.       COMMON /NEXAMP  /NEXAMP
  381.       PI  = DACOS(-1.0D0)
  382. C
  383. C......................................................................
  384. C
  385.       IF (NEXAMP .eq. 1) GOTO 1
  386.       IF (NEXAMP .eq. 2) GOTO 2
  387.       IF (NEXAMP .eq. 3) GOTO 3
  388.       IF (NEXAMP .eq. 4) GOTO 4
  389.       IF (NEXAMP .eq. 5) GOTO 5
  390.       IF (NEXAMP .eq. 6) GOTO 6
  391. C
  392. C......................................................................
  393. C
  394. C     NON OSCILLATING SOLUTION
  395. C
  396.     1 u = dsin(pi*tx)*dsin(pi*ty)*dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
  397.       GOTO 10
  398. C
  399. C......................................................................
  400. C
  401. C     HIGHLY OSCILLATING SOLUTION
  402. C
  403.     2 u = dsin(20.0d0*pi*tx)*dsin(20.0d0*pi*ty)
  404.      +  * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
  405.       GOTO 10
  406. C
  407. C......................................................................
  408. C
  409. C     SOLUTION INCREASING SLOWLY IN TIME
  410. C
  411.     3 u = dexp(tx+ty) * dexp(dsqrt(2.0d0)*t) * 100000.0d0
  412.       GOTO 10
  413. C
  414. C......................................................................
  415. C
  416. C     SOLUTION INCREASING RAPIDLY IN TIME
  417. C
  418.     4 u = dexp(3.0d0*(tx+ty)) * dexp(3.0d0*dsqrt(2.0d0)*t)
  419.       GOTO 10
  420. C
  421. C......................................................................
  422. C
  423. C     POLYNOMIAL IN X AND Y AS SOLUTION
  424. C
  425.     5 u = tx**2 + ty**2 + 2.0d0*t**2
  426.       GOTO 10
  427. C
  428. C......................................................................
  429. C
  430. C     ZERO SOLUTION
  431. C
  432.     6 u = 0.0D0
  433.       GOTO 10
  434. C
  435. C......................................................................
  436. C
  437.    10 CONTINUE
  438.       END
  439. C
  440. C
  441. C****************************************************************
  442. C****************************************************************
  443. C****************************************************************
  444. C****************************************************************
  445. C****************************************************************
  446. C
  447. C
  448.       subroutine F1(u,tx,ty,np)
  449. C
  450. C**********************************************************************
  451. C
  452. C     INITIAL CONDITION U(X,Y,0) = F1(X,Y)
  453. C
  454. C**********************************************************************
  455. C
  456.       integer          np
  457.       DOUBLE PRECISION u(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
  458.       DOUBLE PRECISION PI
  459.       INTEGER          NEXAMP
  460.       COMMON /NEXAMP  /NEXAMP
  461.       PI  = DACOS(-1.0D0)
  462. C
  463. C......................................................................
  464. C
  465.       IF (NEXAMP .eq. 1) GOTO 1
  466.       IF (NEXAMP .eq. 2) GOTO 2
  467.       IF (NEXAMP .eq. 3) GOTO 3
  468.       IF (NEXAMP .eq. 4) GOTO 4
  469.       IF (NEXAMP .eq. 5) GOTO 5
  470.       IF (NEXAMP .eq. 6) GOTO 6
  471. C
  472. C......................................................................
  473. C
  474. C     NON OSCILLATING SOLUTION
  475. C
  476.     1 u = DSIN(PI*tx) * DSIN(PI*ty) * 1000.0D0
  477.       GOTO 10
  478. C
  479. C......................................................................
  480. C
  481. C     HIGHLY OSCILLATING SOLUTION
  482. C
  483.     2 u = DSIN(20.0d0*PI*tx) * DSIN(20.0d0*PI*ty)
  484.       GOTO 10
  485. C
  486. C......................................................................
  487. C
  488. C     SOLUTION INCREASING SLOWLY IN TIME
  489. C
  490.     3 u = Dexp(tx+ty) * 100000.0d0
  491.       GOTO 10
  492. C
  493. C......................................................................
  494. C
  495. C     SOLUTION INCREASING RAPIDLY IN TIME
  496. C
  497.     4 u = Dexp(3.0d0*(tx+ty)) 
  498.       GOTO 10
  499. C
  500. C......................................................................
  501. C
  502. C     POLYNOMIAL IN X AND Y AS SOLUTION
  503. C
  504.     5 u = tx**2 + ty**2
  505.       GOTO 10
  506. C
  507. C......................................................................
  508. C
  509. C     ZERO SOLUTION
  510. C
  511.     6 u = 0.0d0 
  512.       GOTO 10
  513. C
  514. C......................................................................
  515. C
  516.    10 CONTINUE
  517.       END
  518. C
  519. C
  520. C****************************************************************
  521. C****************************************************************
  522. C****************************************************************
  523. C****************************************************************
  524. C****************************************************************
  525. C
  526. C
  527.       subroutine F2(u,tx,ty,np)
  528. C
  529. C**********************************************************************
  530. C
  531. C     INITIAL CONDITION DU/DT(X,Y,0) = F2(X,Y)
  532. C
  533. C**********************************************************************
  534. C
  535.       integer          np
  536.       DOUBLE PRECISION u(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
  537.       DOUBLE PRECISION PI
  538.       INTEGER          NEXAMP
  539.       COMMON /NEXAMP  /NEXAMP
  540.       PI  = DACOS(-1.0D0)
  541. C
  542. C......................................................................
  543. C
  544.       IF (NEXAMP .eq. 1) GOTO 1
  545.       IF (NEXAMP .eq. 2) GOTO 2
  546.       IF (NEXAMP .eq. 3) GOTO 3
  547.       IF (NEXAMP .eq. 4) GOTO 4
  548.       IF (NEXAMP .eq. 5) GOTO 5
  549.       IF (NEXAMP .eq. 6) GOTO 6
  550. C
  551. C......................................................................
  552. C
  553. C     NON OSCILLATING SOLUTION
  554. C
  555.     1 u = 0.0d0
  556.       GOTO 10
  557. C
  558. C......................................................................
  559. C
  560. C     HIGHLY OSCILLATING SOLUTION
  561. C
  562.     2 u = 0.0d0
  563.       GOTO 10
  564. C
  565. C......................................................................
  566. C
  567. C     SOLUTION INCREASING SLOWLY IN TIME
  568. C
  569.     3 u = dsqrt(2.0d0) * dexp(tx+ty) * 100000.0d0
  570.       GOTO 10
  571. C
  572. C......................................................................
  573. C
  574. C     SOLUTION INCREASING RAPIDLY IN TIME
  575. C
  576.     4 u = 3.0d0 * dsqrt(2.0d0) * dexp(3.0d0*(tx+ty))
  577.       GOTO 10
  578. C
  579. C......................................................................
  580. C
  581. C     POLYNOMIAL IN X AND Y AS SOLUTION
  582. C
  583.     5 u = 0.0d0
  584.       GOTO 10
  585. C
  586. C......................................................................
  587. C
  588. C     ZERO SOLUTION
  589. C
  590.     6 u = 0.0d0
  591.       GOTO 10
  592. C
  593. C......................................................................
  594. C
  595.   10  CONTINUE
  596.       END
  597. C
  598. C
  599. C****************************************************************
  600. C****************************************************************
  601. C****************************************************************
  602. C****************************************************************
  603. C****************************************************************
  604. C
  605. C
  606.       subroutine G(u,tx,ty,np,t)
  607. C
  608. C**********************************************************************
  609. C
  610. C     BOUNDARY CONDITION U(X,Y,T) = G(X,Y,T)
  611. C
  612. C**********************************************************************
  613. C
  614.       integer          np
  615.       DOUBLE PRECISION u(0:np,0:np), tx(0:np,0:np), ty(0:np,0:np)
  616.       DOUBLE PRECISION t
  617.       DOUBLE PRECISION PI
  618.       INTEGER          NEXAMP
  619.       COMMON /NEXAMP  /NEXAMP
  620.       PI  = DACOS(-1.0D0)
  621. C
  622. C......................................................................
  623. C
  624.       IF (NEXAMP .eq. 1) GOTO 1
  625.       IF (NEXAMP .eq. 2) GOTO 2
  626.       IF (NEXAMP .eq. 3) GOTO 3
  627.       IF (NEXAMP .eq. 4) GOTO 4
  628.       IF (NEXAMP .eq. 5) GOTO 5
  629.       IF (NEXAMP .eq. 6) GOTO 6
  630. C
  631. C......................................................................
  632. C
  633. C     NON OSCILLATING SOLUTION
  634. C
  635.     1 u(0,0:np) = dsin(pi*tx(0,0:np)) * dsin(pi*ty(0,0:np))
  636.      +          * dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
  637.       u(np,0:np) = dsin(pi*tx(np,0:np)) * dsin(pi*ty(np,0:np))
  638.      +          * dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
  639.       u(1:np-1,0) = dsin(pi*tx(1:np-1,0)) * dsin(pi*ty(1:np-1,0))
  640.      +          * dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
  641.       u(1:np-1,np) = dsin(pi*tx(1:np-1,np)) * dsin(pi*ty(1:np-1,np))
  642.      +          * dcos(dsqrt(2.0d0)*pi*t) * 1000.0d0
  643.       GOTO 10
  644. C
  645. C......................................................................
  646. C
  647. C     HIGHLY OSCILLATING SOLUTION
  648. C
  649.     2 u(0,0:np) = dsin(20.0d0*pi*tx(0,0:np)) 
  650.      +          * dsin(20.0d0*pi*ty(0,0:np))
  651.      +          * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
  652.       u(np,0:np) = dsin(20.0d0*pi*tx(np,0:np)) 
  653.      +          * dsin(20.0d0*pi*ty(np,0:np))
  654.      +          * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
  655.       u(1:np-1,0) = dsin(20.0d0*pi*tx(1:np-1,0))
  656.      +          * dsin(20.0d0*pi*ty(1:np-1,0))
  657.      +          * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
  658.       u(1:np-1,np) = dsin(20.0d0*pi*tx(1:np-1,np)) 
  659.      +          * dsin(20.0d0*pi*ty(1:np-1,np))
  660.      +          * dcos(20.0d0*dsqrt(2.0d0)*pi*t)
  661.       GOTO 10
  662. C
  663. C......................................................................
  664. C
  665. C     SOLUTION INCREASING SLOWLY IN TIME
  666. C
  667.     3 u(0,0:np) = dexp(tx(0,0:np)+ty(0,0:np)) 
  668.      +          * dexp(dsqrt(2.0d0)*t) * 100000.0d0
  669.       u(np,0:np) = dexp(tx(np,0:np)+ty(np,0:np)) 
  670.      +          * dexp(dsqrt(2.0d0)*t) * 100000.0d0
  671.       u(1:np-1,0) = dexp(tx(1:np-1,0)+ty(1:np-1,0)) 
  672.      +          * dexp(dsqrt(2.0d0)*t) * 100000.0d0
  673.       u(1:np-1,np) = dexp(tx(1:np-1,np)+ty(1:np-1,np)) 
  674.      +          * dexp(dsqrt(2.0d0)*t) * 100000.0d0
  675.       GOTO 10
  676. C
  677. C......................................................................
  678. C
  679. C     SOLUTION INCREASING RAPIDLY IN TIME
  680. C
  681.     4 u(0,0:np) = dexp(3.0d0* (tx(0,0:np)+ty(0,0:np)) ) 
  682.      +          * dexp(3.0d0*dsqrt(2.0d0)*t)
  683.       u(np,0:np) = dexp(3.0d0* (tx(np,0:np)+ty(np,0:np)) ) 
  684.      +          * dexp(3.0d0*dsqrt(2.0d0)*t) 
  685.       u(1:np-1,0) = dexp(3.0d0* (tx(1:np-1,0)+ty(1:np-1,0)) ) 
  686.      +          * dexp(3.0d0*dsqrt(2.0d0)*t) 
  687.       u(1:np-1,np) = dexp(3.0d0* (tx(1:np-1,np)+ty(1:np-1,np)) ) 
  688.      +          * dexp(3.0d0*dsqrt(2.0d0)*t) 
  689.       GOTO 10
  690. C
  691. C......................................................................
  692. C
  693. C     POLYNOMIAL IN X AND Y AS SOLUTION
  694. C
  695.     5 u(0,0:np) = tx(0,0:np)**2 + ty(0,0:np)**2 + 2.0d0*t**2
  696.       u(np,0:np) = tx(np,0:np)**2 + ty(np,0:np)**2 + 2.0d0*t**2
  697.       u(1:np-1,0) = tx(1:np-1,0)**2 + ty(1:np-1,0)**2 + 2.0d0*t**2
  698.       u(1:np-1,np) = tx(1:np-1,np)**2 + ty(1:np-1,np)**2 + 2.0d0*t**2
  699.       GOTO 10
  700. C
  701. C......................................................................
  702. C
  703. C     ZERO SOLUTION
  704. C
  705.     6 u(0,0:np) = 0.0d0
  706.       u(np,0:np) = 0.0d0
  707.       u(1:np-1,0) = 0.0d0
  708.       u(1:np-1,np) = 0.0d0
  709. C
  710. C......................................................................
  711. C
  712.    10 CONTINUE
  713.       END
  714.