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