# Rcjp's Weblog

## February 14, 1991

### BSc Fortran Numerics Course – Solving Equations

Filed under: physics — rcjp @ 11:06 am

We did some very basic differential equation solving, simultaneous equations and curve fitting…

``````
CCCC  ******************************************************************
CCCC  **        SOLUTION OF FIRST ORDER DIFFERENTIAL EQUATION         **
CCCC  **                USING EULERS METHODS                          **
CCCC  ******************************************************************

EXTERNAL FUNC

REAL XO,XOUT,Y,H,FUNC
DATA XO,XOUT,Y/0.0,3.0,1.0/

PRINT *,'ENTER STEP SIZE '

CALL EULER (XO,XOUT,H,Y,FUNC)

PRINT *,'Y VALUE AT X=3.0 IS ',Y
END

REAL FUNCTION FUNC (X,Y)
REAL X,Y
FUNC=-Y
END

SUBROUTINE EULER (XO,XOUT,H,Y,FUNC)
REAL I,XO,XOUT,H,Y
DO 100 I=XO,XOUT,H
Y=Y+H*FUNC (XO,Y)
100 CONTINUE
END
``````

and

``````
CCCC  ******************************************************************
CCCC  **        SOLUTION OF FIRST ORDER DIFFERENTIAL EQUATION         **
CCCC  **                USING RUNGE-KUTTA METHODS                     **
CCCC  ******************************************************************

EXTERNAL FUNC

REAL XO,XOUT,Y,H,FUNC
DATA XO,XOUT,Y/0.0,3.0,1.0/

PRINT *,'ENTER STEP SIZE '

CALL RK4 (XO,XOUT,H,Y,FUNC)

PRINT *,'Y VALUE AT X=3.0 IS ',Y

END

REAL FUNCTION FUNC (X,Y)
REAL X,Y
FUNC=-Y
END

SUBROUTINE RK4 (XO,XOUT,H,Y,FUNC)
REAL X,XO,XOUT,H,Y,C1,C2,C3,C4
DO 100 X=XO,XOUT,H
C1=H*FUNC (X,Y)
C2=H*FUNC (X+0.5*H,Y+0.5*C1)
C3=H*FUNC (X+0.5*H,Y+0.5*C2)
C4=H*FUNC (X+H,Y+C3)
Y=Y+(C1+2*C2+2*C3+C4)/6
100 CONTINUE
END
``````

## Simultaneous Equations

``````
CCCC  ******************************************************************
CCCC  **            SOLUTION OF SIMULTANEOUS EQUATIONS                **
CCCC  **              BY GAUSSIAN ELIMINATION METHOD                  **
CCCC  **            WITH PARTIAL PIVOTAL CONDENSATION                 **
CCCC  **                AND CHECK FOR SINGULARITY                     **
CCCC  ******************************************************************

CCCC  Initialise variables
INTEGER I,J,K,N,L
REAL A,B,X,M,F,TOL,C,D,S
DIMENSION A(20,20),X(20),B(20)
COMMON A,B,X

CCCC  Enter equations

PRINT*,'ENTER NO. OF EQUATIONS AND TOLERANCE '
DO 10 I=1,N
PRINT114,'ENTER COEFFICIENTS OF EQU.',I
10 CONTINUE
PRINT113,'EQUATIONS ENTERED'
CALL PR(N)

CCCC  PARTIAL PIVOTAL CONDENSATION

DO 5 J=1,N
L=J
DO 2 I=1,N
IF (A(I,J).GT.A(L,J)) L=I
2   CONTINUE
DO 3 K=1,N
C=A(J,K)
A(J,K)=A(L,K)
A(L,K)=C
3   CONTINUE
D=B(J)
B(J)=B(L)
B(L)=D
5 CONTINUE
PRINT113,'PARTIAL PIVOTAL CONDENSATION'
CALL PR(N)

CCCC  Gaussian elimination method

DO 50 K=1,N-1
DO 40 I=K+1,N
IF (ABS(A(K,K)).LT.TOL) STOP
M=A(I,K)/A(K,K)
S=0
DO 30 J=K,N
A(I,J)=A(I,J)-M*A(K,J)
S=S+A(I,J)**2
30     CONTINUE
B(I)=B(I)-M*B(K)
IF (S.LT.TOL) THEN
PRINT113,'SINGULARITY ERROR'
STOP
ENDIF
40   CONTINUE
50 CONTINUE
PRINT113,'TRIANGULAR FORM BEFORE BACK SUBSTITUTION'
CALL PR(N)

CCCC  Back substitution & output of results

PRINT113,'SOLUTIONS'
X(N)=B(N)/A(N,N)
PRINT115,N,X(N)
DO 110 J=N-1,1,-1
DO 100 L=J+1,N
F=F+A(J,L)*X(L)
100   CONTINUE
X(J)=1/A(J,J)*(B(J)-F)
PRINT115,J,X(J)
F=0
110 CONTINUE
113 FORMAT (/,1X,A)
114 FORMAT (1X,A,I1,6X)
115 FORMAT (1X,'X(',I1,') = ',F9.5)
END

CCCC  Display equations

SUBROUTINE PR(N)
INTEGER I,J,N
COMMON A(20,20),B(20),X(20)
DO 210 I=1,N
DO 200 J=1,N
PRINT215,A(I,J)
200   CONTINUE
PRINT216,B(I)
210 CONTINUE
215 FORMAT (1X,F9.4,4X,\)
216 FORMAT (1X,'=',2X,F9.4)
END
``````

We were obviously using the Salford FTN77 compiler since there is a log file

``````
SALFORD UNIVERSITY FTN77/386 VER. 1.67      A:\GAUSSPC.FOR    15:29:42    Tuesday, 22 January  1991
COMPILER OPTIONS:  LISTING ANSI INTL NODCLVAR NOMAP NOCHECK LOGS DYNM OFFSET LGO NOPAGETHROW NOSILENT
NO_OPTIMISE NOWEITEK
``````

I can remember having to go to the machine room in the Newton building to get the listings from the old 132 column printer (everywhere else was just epson 850’s which couldn’t print that wide.)

## Finite Difference

``````
CCCC  ******************************************************************
CCCC  **   CALCULATES THE TEMPERATURE IN A BLOCK OF HEATED MATERIAL   **
CCCC  **   USES INITIAL GUESSES AND RELAXATION TECHNIQUE APPLIED TO   **
CCCC  **   THE FINITE DIFFERENCE EQUATIONS AND DRAWS A GRAPH OF THE   **
CCCC  **   NUMBER OF ITERATIONS AGAINST THE RELAXATION COEFF. ALPHA   **
CCCC  ******************************************************************

CCCC  Block is heated ...
CCCC      500 500 500 500
CCCC      100 ??? ??? 100
CCCC      100 ??? ??? 100
CCCC      100 100 100 100

REAL T,RES,RESMX,DECT,W,WMAX,WMIN,R,TT
INTEGER I,J,M,N,TMMIN,TMMAX,TNMIN,TNMAX,IT,PTS,NUM
DIMENSION T(10,10),TT(10,10),RES(10,10),PTS(50),W(50)

DATA  ((TT(I,J),I=1,4),J=1,4),TMMIN,TMMAX,TNMIN,TNMAX/
1  4*500.0,100.0,2*300.0,2*100.0,2*200.0,5*100.0,1,4,1,4/

PRINT *,'ENTER RELAXATION COEFFICIENT BOUNDARIES '

DO 250 R=WMIN,WMAX,(WMAX-WMIN)/40.0

DO 5 I=1,4
DO 4 J=1,4
T(I,J)=TT(I,J)
4 CONTINUE
5 CONTINUE

DO 200 IT=0,200

CCCC  FIND RESIDUALS

RESMX=0.0
DO 20 M=TMMIN+1,TMMAX-1
DO 10 N=TNMIN+1,TNMAX-1
RES (M,N)=T(M+1,N)+T(M-1,N)+T(M,N+1)+T(M,N-1)-4.0*T(M,N)
IF (ABS(RES(M,N)).GT.RESMX) RESMX=RES(M,N)
10     CONTINUE
20   CONTINUE

CCCC  ALTER TEMPERATURE OF ALL POINTS BY RESMX/4.0*R

DECT=RESMX/4.0*R
DO 40 M=TMMIN+1,TMMAX-1
DO 30 N=TNMIN+1,TNMAX-1
T(M,N)=T(M,N)+DECT
30     CONTINUE
40   CONTINUE
IF (ABS(RESMX).LT.0.1 .AND. IT.GT.0) GOTO 230
200 CONTINUE
C  230 PRINT *,'REL COEFF',R,'  STEPS ',IT
230 NUM=NUM+1
W(NUM)=R
PTS(NUM)=IT
250 CONTINUE
300 FORMAT (1X,I3,1X,12(F6.1))

CALL SCRSET
CALL PPOINT (NUM,W,PTS)

END

CCCC  Initialise screen and clear it.

SUBROUTINE SCRSET
INCLUDE 'F77PC'
CALL INITSCREEN
CALL SETSCREENMODE (EGAGRAPH)
CALL CLRVIDEO
END

SUBROUTINE PPOINT (NPTS,X,Y)
INCLUDE 'F77PC'
INTEGER NPTS,Y,N,I,XP,YP
REAL X,XMAX,YMAX,SCALEX,SCALEY,XMIN,YMIN
DIMENSION Y(50),X(50)
CCCC  Draw axis

CALL DRAW (20,330,620,330,WHITE)
CALL DRAW (20,330,20,10,WHITE)

CCCC  Mark points
XMIN=1000
YMIN=1000
DO 390 N=1,NPTS
IF (X(N).GT.XMAX) XMAX=X(N)
IF (X(N).LT.XMIN) XMIN=X(N)
IF (Y(N).GT.YMAX) YMAX=Y(N)
IF (Y(N).LT.YMIN) YMIN=Y(N)
390 CONTINUE

SCALEX=600.0/(XMAX-XMIN)
SCALEY=320.0/(YMAX-YMIN)
DO 400,I=1,NPTS
XP=INT((X(I)-XMIN)*SCALEX+20)
YP=INT(330-(Y(I)-YMIN)*SCALEY)-2
CALL DRAW (XP-2,YP,XP+2,YP,WHITE)
CALL DRAW (XP,YP-2,XP,YP+2,WHITE)
400 CONTINUE
END
``````

## Least Squares Fitting

``````
CCCC  ******************************************************************
CCCC  **      DRAWS A GRAPH OF DATA HELD IN FILE 'GDATA.DAT'          **
CCCC  **               AND DRAW LEAST SQUARES FIT                     **
CCCC  ******************************************************************

CCCC  Data file includes the number of points
CCCC  and X,Y coordinates

INCLUDE 'F77PC'
INTEGER NPTS,I
DIMENSION X(30),Y(30)

OPEN (5,FILE='A:\LSQ\GDATA.DAT',STATUS='OLD')
DO 10 I=1,NPTS
10 CONTINUE

CALL SCRSET
PRINT *,'INTERCEPT=',INTCPT
CLOSE (5)
END

CCCC  Calculation of least squares fit to data
CCCC  with coordinates held in array X,Y.
CCCC  on the y -axis

INTEGER NPTS,I
DIMENSION X(*),Y(*)
DO 250 I=1,NPTS
SX=SX+X(I)
SXX=SXX+X(I)*X(I)
SXY=SXY+X(I)*Y(I)
SY=SY+Y(I)
250 CONTINUE
DEL=NPTS*SXX-(SX*SX)
INTCPT=(SXX*SY-SXY*SX)/DEL
END

CCCC  Initialise screen and clear it.

SUBROUTINE SCRSET
INCLUDE 'F77PC'
CALL INITSCREEN
CALL SETSCREENMODE (EGAGRAPH)
CALL CLRVIDEO
END

CCCC  Plots the data, draws the axis
CCCC  and the line of best fit

INCLUDE 'F77PC'
INTEGER I,NPTS,XP,YP,Y1,Y2
DIMENSION X(*),Y(*)

CCCC  Draw axis
CALL DRAW (20,330,620,330,WHITE)
CALL DRAW (20,330,20,10,WHITE)

CCCC  Mark points
SCALEX=600/X(NPTS)
SCALEY=320/Y(NPTS)
DO 300,I=1,NPTS
XP=INT(X(I)*SCALEX+20)
YP=INT(330-Y(I)*SCALEY)
CALL DRAW (XP-2,YP,XP+2,YP,WHITE)
CALL DRAW (XP,YP-2,XP,YP+2,WHITE)
300 CONTINUE

CCCC  Draw best fit line
Y1=INT(330-INTCPT*SCALEY)