Get documents about "Ke-42.4520 Process Modeling - methods and tools FORTRAN hands on
Document Sample
Ke-42.4520 Process Modeling -
methods and tools
FORTRAN hands on
D.Sc.(Tech.) Kaj Jakobsson
Helsinki University of Technology
Chemical Engineering and Plant Design
POB 6100
FIN-02015 HUT
FINLAND
email: Kaj.Jakobsson@hut.fi
Fortran history
FORTRAN I
Development 1954-57 by an IBM team
FORTRAN II, III, IV and FORTRAN 66
Development during 1958-1966
FORTRAN 77
Formally “outdated” years ago, compilers for
FORTRAN 77 are still used today
Fortran 90/95
Fortran 90
Adds many extensions to FORTRAN 77.
The language in its present form is
competitive with computer languages
created later (e.g. C).
Fortran 95
Added some minor improvements to the
Fortran 90 standard.
FORTRAN ADVANTAGES
a) Scientifically oriented
b) Better optimized code
c) A lot of existing code
d) Easy to learn
e) Efficient mathematics
f) Easy to use and robust
g) Good diagnostics
Our Standard (HUT Chem. Eng)
Strict FORTRAN 77 with some widely
accepted extensions
DO …END DO
NAMELIST directed input
INCLUDE
Motivation: Maximum portability
Program layout
Column 1 C or * to begin a comment line
Otherwise Columns 7 to 72 are FORTRAN
stament, or part of one
Columns 1..5 may continue a label unless 6
has continuation symbol.
Columns 73 and onwards are comments
Layout examples
C123456
C This is a comment
C Below a continued line
A = ALFA + BETA +
& GAMMA
(Blanks are not important in FORTRAN)
Below an example of a label
1000 CONTINUE
Layout pitfall
You write a variable so that it extends over 73
column. The variable name changes and
messes up you calculations
column 73
stage
Value names
In FORTRAN 77 the name must start with a
letter.
Names can have up to 6 symbols (letter or
symbols)
(New standards are more free in this
respect)
Data types
The basic data types of FORTRAN 77 are:
integers, floating-point numbers of two
kinds (REAL and DOUBLE PRECISION),
logicals and character strings.
Implicit statement
IMPLICIT DOUBLE PRECISION (A - H, O - Z)
This means that all values named so that they are
starting with A-H or O-Z are double precision, I-
N are integers.
Often used in old codes.
DOUBLE PRECISION
We use in practice only double precision for
floating point calculations
C23456
DOUBLE PRECSION MODELV
If vector or array
DIMENSION MODELV(100)
Or
DIMENSION MODELV(100,100)
Etc…
DOUBLE PRECISION
Double precision real constants
2.0 2.71824 -1 .1D-1
-0.1D-3 6.025D23
Some aritmethic
FORTRAN is made for computing
Operations has same priority as
mathematics, operations with same priority
are calculated from left to right
Some arithmetic
Addition A = 3.0 + 4.0
Subtraction B = 15D2 – 12
Multiplying C=A*B
Division D=A/B
Exponent E = X**Y
Bracets ( ) are used as in mathematics
Note In FORTRAN you cannot use more than one
operator between two values for example 9.00 +
(-3)
INTEGER
C23456
INTEGER MODELV
If vector or array
DIMENSION MODELV(100)
Or
DIMENSION MODELV(100,100)
Etc…
INTEGER
Integer constants
4 -100 +1024 35000
CHARACTER
C23456
CHARACTER EMPTY * 20
CHARACTER PHASE(MAXDFL) * 20
…
EMPTY = „The string'
LOGICAL
C23456
LOGICAL IHIT
Can be .TRUE. or .FALSE.
Relational Operators F77
.lt. Less than
.le. Less or equal
.eq. Equal
.ne. Not equal
.gt. Greater than
.ge. Greater or equal
The periods are obligatory
Logical expressions
.NOT. Logical expression
If logical expression .TRUE. result .FALSE. and
vice versa
Logical expression .AND. Logical expression
The result is true if both are .TRUE. else .FALSE.
Logical expression .OR. Logical expression
The result is true if either or both of the logical
expressions are .TRUE.
Priority of logical expressions
()
**
Or /
+ or –
// Character operation
.EQ. .NE. .GT. .LT. .GE. .LE.
.NOT. .AND.
.OR.
Logical IF STATEMENT
IF(I.EQ.1) VARIAB=4.5
BLOCK IF STATEMENT
IF(I.EQ.0) THEN
A=4.5
B=TAN(A)
END IF
IF WITH ALTERNATIVES
IF(I.EQ.1) THEN
B =1.0
ELSE IF(I.EQ.2) THEN
B=0.0
ELSE IF(LOGICAL EXPRESSION)
ELSE IF(LOGICAL EXPRESSION)
…..
ELSE
B=-1.0D25
END IF
DO LOOP
DO 10 I=1,100
ARR(I)=5.0*ARR(I)
10 CONTINUE
DO loops are used to repeat the same set of
instructions many times over, for example,
when accessing all the elements of an
array one at a time.
DO LOOP
TOT=0.0
DO 20 I=1,100,2
TOT=TOT+ARR(I)**3.0
20 CONTINUE
This loop adds the cubes of every second
element in ARR to TOT, which is set to
zero initially.
DO LOOP
CALL WRC1('BMTRCL Warning -- Unknown method for KL-values')
DO 110 I=1,NC-1
DO 100 J=I+1,NC
KL(I,J)=1.0D-5
KL(J,I)=1.0D-5
100 CONTINUE
KL(I,I)=0.D0
110 CONTINUE
KL(NC,NC)=0.D0
Do loops can be inside another loop also
Test after (Repeat)
F77 have not an explicit implementation for
this structure. However it can be done
label Labelled statement
Some statements
IF(.NOT. Some condition) GO TO label
Test before (WHILE)
start label IF(.NOT. Some condition) GO TO final label
Some statements
GO TO start label
Final label final statement
SUBROUTINES
It is easier to gather algorithms into
independent modules which can be
shared through program libraries
These algorithms can be mathematical,
physical properties calculations, unit
operations, data sorting etc…
SUBROUTINES
In FORTRAN 77 we have FUNCTIONs and
SUBROUTINES
These modules can be added into programs
without changes
They share data with other program units
only with parameter lists and COMMON
areas
USER DEFIEND FUNCTIONS
Every function must have a type and the function name must be
assigned a value within the function itself CROOT. The arguments are
past by a parameter list (X,Y)
DOULE PRECISION FUNCTION CROOT(X,Y)
C ***************************************************
DOUBLE PRECISION X, Y
IF((X*Y).GT.0) THEN
CROOT=(X*Y)**(1/3.0)
ELSE
CROOT=0.0D0
END IF
RETURN
END
The function is called by actual arguments
…
Z = CROOT(A,B)
…
USER DEFIEND FUNCTIONS
Functions can be also
INTEGER FUNCTIONS …
REAL FUNCTIONS …
LOGICAL FUNCTIONS …
LOGICAL FUNCTION LOTEST(X, T)
C
C Test if X is integer * T
C
IMPLICIT DOUBLE PRECISION (A - H, O - Z)
C
A=X/T
IA = INT(A)
C
IF ((A - DBLE(IA)) .LT. 1.0D-7) THEN
LOTEST = .TRUE.
ELSE
LOTEST = .FALSE.
END IF
C
RETURN
C
END
SUBROUTINE
SUBROUTINE NAME (DUMMY ARGUMENTS)
LOTS OF STATEMENTS
perhaps some RETURNs
END
example
C23456
SUBROUTINE ADDUP(BARRAY,TOTAL)
C ***********************************************
DOUBLE PRECISION BARRAY(100), TOTAL
C
DO 10 I=1,100
TOTAL=TOTAL+BARRAY(I)
10 CONTINUE
C THE FINAL TWO LINES OF A SUBROUTINE MUST BE
C"RETURN" AND "END“ RESPECTIVELY.
RETURN
END
The subroutine is called
CALL ADDUP(B,T)
Example adjustable size array
C23456
SUBROUTINE ADDUP(N, BARRAY,TOTAL)
C ***********************************************
DOUBLEPRECISION BARRAY(*), TOTAL
INTEGER N
DO 10 I=1,N
TOTAL=TOTAL+BARRAY(I)
10 CONTINUE
C THE FINAL TWO LINES OF A SUBROUTINE MUST BE
C"RETURN" AND "END“ RESPECTIVELY.
RETURN
END
Main program
PROGRAM TEST
DOUBLEPRECISION ARRAY, TOTAL
INTEGER NTOT
DIMENSION ARRAY(1000)
NTOT=100
DO 10 I=1, NTOT
ARRAY(I) = 1.0D0
10 CONTINUE
CALL ADDUP(NTOT, ARRAY, TOTAL)
PRINT*, 'TOTAL:', TOTAL
PAUSE
END
Visualisation of the computer
memory
PROGRAM TEST
DOUBLEPRECISION ARRAY, TOTAL
INTEGER NTOT 1 1000
DIMENSION ARRAY(1000)
NTOT=100
DO 10 I=1, NTOT
ARRAY(I) = 1.0D0
10 CONTINUE
CALL ADDUP(NTOT, ARRAY, TOTAL)
PRINT*, 'TOTAL:', TOTAL
PAUSE
END
Visualisation of the computer
memory
1 1000
C23456
SUBROUTINE ADDUP(N, BARRAY,TOTAL)
C ***********************************************
DOUBLEPRECISION BARRAY(*), TOTAL
INTEGER N
DO 10 I=1,N
TOTAL=TOTAL+BARRAY(I)
10 CONTINUE
C THE FINAL TWO LINES OF A SUBROUTINE MUST BE
C"RETURN" AND "END“ RESPECTIVELY.
RETURN
END
Visualisation of the computer
memory
1 1000
It is easy mess up in the calls with wrong
types (for example the integer and double
precision variables are of different length)
or wrong array size dimensioning.
Arrays
DOUBLEPRECISION A, B
DIMENSION A(3, 3), B(5, 5, 5)
To make flexible FORTRAN 77 programs
arrays are first tricky to understand.
Arrays
Useful idea when you declare an array in
your program: you declare an working
space with an certain length.
DOUBLEPRECISION A
DIMENSION A(3, 3)
1,1 2,1 3,1 1,2 2,2 2,3 1,3 2,3 3,3
In Fortran the first index runs fastest
Why TOTAL 3? Should be 4
PROGRAM TEST
DOUBLEPRECISION ARRAY, TOTAL C23456
INTEGER NTOT SUBROUTINE ADDUP(N, BARRAY,TOTAL)
DIMENSION ARRAY(3,3) C ***********************************************
DOUBLEPRECISION BARRAY(N,N), TOTAL
NTOT=2 INTEGER N
C
DO 10 I=1, NTOT TOTAL = 0.0D0
DO 110 J=1,NTOT DO 10 I=1,N
ARRAY(I,J) = 1.0D0 DO 110 J=1,N
110 CONTINUE TOTAL=TOTAL+BARRAY(I,J)
10 CONTINUE 110 CONTINUE
10 CONTINUE
CALL ADDUP(NTOT, ARRAY, TOTAL) C THE FINAL TWO LINES OF A SUBROUTINE MUST BE
C"RETURN" AND "END“ RESPECTIVELY.
PRINT*, 'TOTAL:', TOTAL RETURN
PAUSE END
END
Answer
DOUBLEPRECISION ARRAY
DIMENSION A(3, 3)
Main program fills these memory locations
1,1 2,1 3,1 1,2 2,2 2,3 1,3 2,3 3,3
Subroutine uses these memory locations
Text books says:
Only last index flexible in subroutines
C23456
SUBROUTINE ADDUP(N, BARRAY,TOTAL)
C ***********************************************
DOUBLEPRECISION BARRAY(3,*), TOTAL
INTEGER N
C
TOTAL = 0.0D0
DO 10 I=1,N
DO 110 J=1,N
TOTAL=TOTAL+BARRAY(I,J)
110 CONTINUE
10 CONTINUE
C THE FINAL TWO LINES OF A SUBROUTINE MUST
BE C"RETURN" AND "END“ RESPECTIVELY.
RETURN
END
HOWEVER
The problem can be solved with defining a
working area in the main program or an
„envelope‟ routine to do these tasks in a
subroutine.
It is useful to have a flexible program. We
need different number of components,
reactors, modelling stages ETC…
Example solution
PROGRAM TEST SUBROUTINE ARR(NTOT, TOTAL, ARRAY)
DOUBLEPRECISION WA DOUBLEPRECISION ARRAY, TOTAL
INTEGER IT, IA, NTOT INTEGER NTOT
DIMENSION WA(20000)
DIMENSION ARRAY(NTOT,NTOT)
C Define starting indexes to working area vector WA
C for the vectors below DO 10 I=1, NTOT
IT = 1 DO 110 J=1,NTOT
IA = 1 + IT ARRAY(I,J) = 1.0D0
NTOT = 5 110 CONTINUE
10 CONTINUE
CALL ARR(NTOT, WA(IT), WA(IA))
CALL ADDUP(NTOT, ARRAY, TOTAL)
PRINT*, 'TOTAL:', WA(IT)
PAUSE
PRINT*, 'TOTAL:', TOTAL
END PAUSE
END
Example solution
SUBROUTINE ADDUP(N, BARRAY,TOTAL)
C ***********************************************
DOUBLEPRECISION BARRAY(N,N), TOTAL
INTEGER N
C
TOTAL = 0.0D0
DO 10 I=1,N
DO 110 J=1,N
TOTAL=TOTAL+BARRAY(I,J)
110 CONTINUE
10 CONTINUE
C THE FINAL TWO LINES OF A SUBROUTINE MUST
BE C"RETURN" AND "END“ RESPECTIVELY.
RETURN
END
EXTERNAL
You can pass the names of the subroutines using the EXTERNAL
statement. This is useful in multipurpose codes like mathematical
solvers.
…
EXTERNAL TEST
…
CALL KOE(TEST,…
This implies that some where exists a called
SUBROUTINE TEST
It has the same parameter list as the subroutine
that is called by KOE
COMMON BLOCKS
With a COMMON block you can transport
data between different modules in a
program
It is best used for transporting constant type
data, like data-base values to other
program units
COMMON BLOCKS
COMMON / HENRY / HENRYA(NOFC), HENRYB(NOFC), HENRYC(NOFC)
By adding this line into your PROGRAMs, FUNCTIONs or
SUBROUTINEs you can share the data all over your
program in the particular program unit.
The particular COMMON must be exactly the same size
everywhere in the program
You can assign data to the variables or use the BLOCK
DATA
BLOCK DATA
BLOCK DATA
……
PARAMETER (NOFC = 40)
COMMON / HENRY / HENRYA(NOFC), HENRYB(NOFC), HENRYC(NOFC)
…
DATA HENRYA / NOFC * 0.0D0 /
DATA HENRYB / NOFC * 0.0D0 /
DATA HENRYC / NOFC * 0.0D0 /
….
END
Simple I/O system input
OPEN a text file for input:
NIN = 20
OPEN(NIN, FILE =„C:\TEMP\WATCOM\KOE.DAT')
READ from this file a NAMELIST -input:
READ(NIN, KOE)
NIN is called “channel number” and must be transported all along the program where
ever needed to connect the input to the file.
In this course the NAMELIST input is strongly recommended
NAMELIST
Common extension among FORTRAN 77 compilers,
incorporated into the Fortran 90 standard.
Provides a convenient alternative to writing one's own
parser for free-form "keyword=value" input,
Its main drawbacks are: - the external representation of
NAMELIST data varies between implementations. Some
require / as a terminator, some want &END or $END.
Some ignore column 1 on input, others don't. On output,
some prefer $ to & as a prefix. - there is no guarantee
that any given FORTRAN 77 compiler will support
NAMELIST.
Example of NAMELIST
&KOE
NVAR = 2
A= 1, 50,
B =1, -0.5,
C=1, -0.5,
D=1, -0.5,
XMIN= 0.0, 0.0,
XMAX= 100.0, 100.0,
XVAR= 1.0, 1.0,
TOL =1.0d-6,
IPR = 3
&END
NAMELIST RULES
Detailed rules for namelist input can be
found elsewhere
At HUT for example from Flowbat –manual
http://www.tkk.fi/Units/ChemEng/WebFlowb
at/FlowManual/Manuaali.doc
Output
NCAL = 21
OPEN(NCAL, FILE = 'C:\TEMP\WATCOM\KOE.OUT')
NCAL is called “channel number” and must be transported all along the program where
ever needed to connect the output to the file.
Output and FORMATs
WRITE(NCAL,100)
WRITE(NCAL,110) NVAR
WRITE(NCAL,*) ' XVAR: '
WRITE(NCAL,120) (XVAR(I),I=1,NVAR)
WRITE(NCAL,130) TOL
WRITE(NCAL,140) RES
WRITE(NCAL,150) LL
100 FORMAT(' ** Variables and parameters before solver routine **')
110 FORMAT(' NVAR:', I5)
120 FORMAT(5G16.4)
130 FORMAT(' TOL:', G16.4)
140 FORMAT(' RES:', G16.4)
150 FORMAT(' LL:', I5)
** Variables and parameters before solver routine **
NVAR: 2
XVAR:
1.000 1.000
XMAX:
100.0 100.0
XMIN:
0.0000E+00 0.0000E+00
TOL: 0.1000E-05
RES: 1000.
LL: 0
Compiling and linking
A compiler translates FORTRAN codes into binary
codes called the object codes. The files
containing these translated codes generally
have the extension of ".o" or ".obj.
A linker program reads the object file(s), identifies
all routines needed, searches for them in the
given object files, collects all pieces of
instruction needed for the computer to carry out
the given task, and combines them into an
executable file
PROGRAMMING RULES
a) Clarity is more important than efficiency
b) Generalize your code
c) Write modular programs
d) Document your program
e) Write standard FORTRAN
f) Improve your layout
g) Program defensively
h) Use appropriate algorithms
How to do in practice?
There exist a huge amount of ready code for
open source mathematical code in the
web
Your company or your community has
probably already available a lot of useful
stuff models, physical property
calculations etc.
How to do in practice?
You combine the already existing code with
your own modelling work.
National Institute of Standards and
Technology
The Transactions on Mathematical Software (TOMS)
http://math.nist.gov/toms/
Decision tree GAMS (Guide to Available Mathematical Software)
http://gams.nist.gov/
Netlib
Netlib is a collection of mathematical
software, papers, and databases.
http://www.netlib.org/
Sourceforge.net
SourceForge.net is the world's largest Open
Source software development web site.
http://sourceforge.net/
Zuse Institute Berlin
The Zuse Institute Berlin is a research
institute for applied mathematics and
computer science.
http://www.zib.de/
Mathematical Tools in Flowbat
Our inhouse process flowsheet simulator Flowbat has a large
library of mathematical tools available.
Ordinary Differential Equation (ODE) –solvers, for many
different cases
Differential Algebraic Equation (DAE) -solver
Non-linear algebraic equation solvers, for many different types
Matrix operations
You can get them to work for you by making a small FORTRAN
program and link it to our library routines.
Some non-linear equation solvers
http://bucky.stanford.edu/numericalmethods/
PUBCODE/chapter_5.htm
There are many more…
Numerical Recipes
Excellent book, commercial
Electronic form with sources
http://www.nr.com/
Press, William H.
Numerical recipes in Fortran 77 : the art of scientific computing - 2nd. ed. New
York (NY), 1996. - 933 s. –
Press, William H.
Numerical recipes : the art of scientific computing New York, Cambridge
University Press 2007
Exercise I
Two polynoms
Y =1 + X + X2 + X3
Y = 50 -0.5 X – 0.5 X2 - 0.5 X3
Exercise I
The course CD contains the file equs10.for, which contains the
Newton –solver
The model is in file user1.for
The template input file is in KOE.DAT
Explore the solver routines, mainly the parameter list. Explore the
input file, complete the model file, build the program and solve the
exercise
Newton‟s method
In the Newton’s method the new value for the vector of variables are xk+1 are
reached in the following way
x k 1 x k S x k 1
x k +1 = J k Fk
-1
F
J k
x k
T
F F1 ...FN
T
T T
xT x1 ...x N
T
T T
The iteration continues until a desired termination criteria has been met.
Typically a vector norm has to be satisfied
F 2
TOL
Or a maximum number of iteration k > kmax
Graphical solution
50
45
40
35
30
25
20
15
10
5
0
-5 0 1 2 3 4 5 6
-10
Exercise II
Solve
dY(1)/dt = 1195.0*Y(1) - 1995.0 * Y(2)
dY(2)/dt = 1197.0*Y(1) - 1997.0 * Y(2)
Initial values: Y(1) = 2.0, Y(2) = -2.0
Integration t 0 -> 0.1
analytical solution dY(1)/dt (0.1)=8.187,
dY(2)/dt (0.1)=4.912
Exercise II
On the course CD you can find the files
Lapack.for
Lsoda.for
Odepack.for
These files makes the ode-solver we use
User1.for is the model file
Explore the documentation of the solver, complete the
model, build the program and solve the exercise
Introduction ODEs
Remind you of the basic ideas of the numerical solving of Ordinary
Differential Equations ODEs.
Remind you of the concepts and definitions to get a grip on the
large variety of methods for solving ODEs.
Where ODEs
Simple batch reactor model
k1
Reaction A B P
1 dnA
k1C ACB
V dt
1 dnB
k1C ACB Time is the independent
V dt
variable!
1 dnP
k1C ACB
V dt
assuming there is no volume change 1 dnI dCi
V dt dt
Where ODEs
Simple model of a plug flow reactor
k1
A B P
dCA
k1C ACB
d
Space time or the length coordinate
dCB is the independent variable!
k1C ACB
d
dCP
k1C ACB
d
where is the space time A dL
d
Q
Where ODEs
Modelling Emergency Relief for Processes at Near Critical Conditions
Juha-Pekka Pokki, Juhani Aittamaa, Kari I. Keskinen, Markku Hurme
Computers and Chemical Engineering (1999) S399-S401
dmi dmPOL C
mi ,in mi ,out r i , i 1..c mPOL,in mPOL,out r
dt dt i 1
i
dmCAT dH
mCAT ,in mCAT ,out Ein Eout
dt dt
dHCAT dH POL
ECAT ,in ECAT ,out EPOL,in EPOL,out
dt dt
The independent variable time!
dUTOT dH dH POL dH CAT
dt dt dt dt Juha-Pekka used the Euler’s method to solve this
system.
Where ODEs
To solve difficult steady state problems ’dynamically’ through transient
states mass transfer stage j
0 V j M tj
V
UV Y M V
j ijV ij
UV H j EV
I V
j j
0
Yij Rij
0 YNC , j S V
I
Liquid Nj Vapour I I
j
0 d
X ij Qij
QLj Ej QVj 0 dt T jI E I
j
0 L j M tjL
UL X M L
j ij ij
H j EL
l L
U
N L
j j
0 ij Rij
0 N nc , j S L
j
Where ODEs
Dynamic simulation of processes
The importance and use of dynamic process simulation
Startup and shut down studies, Effects of process parameters, Hazard analysis,
Strategies for multi-product plants, Response of equipment failure, Operation
training, Operability Studies, Control strategies.
dy
The general simulation f ( y, z, t , P)
equations dt
0 g ( y, z, t, P)
Initial conditions must be y (a ) y 0
satisfied
Sometimes the exact corresponding values z(a) that satisfy the equations must be known
Numerical solution of ODE –
initial value problems
Problems involving ordinary differential equations (ODEs) can always be reduced to
study of sets of first order differential equations.
dy
f ( x, y ), y(a ) y 0 ; x (a, b); y R n
dx
How to reduce, an example: d2y dy
q( x ) r ( x )
dx 2 dx
dy
This equation can be rewritten as two first order z (x )
equations: dx
dz
r ( x ) q( x ) z ( x )
dx
The usual choice to introduce new variables is to let them be just derivatives of each other but let
common sense be your guide.
Some definitions
The underlying idea of any routine for solving the initial value problem is always:
Rewrite dy and dx as finite steps y and x and and multiply the equations by x.
This gives algebraic formulas for the change in functions when the independent
variable is changed by x.
The method is said to be explicit if the algebraic formula is explicit
in yn
The Euler’s method is yn 1 yn x f ( yn )
explicit:
The method is said to be implicit if the algebraic formula is implicit
in yn
The backward Euler is yn 1 yn x f ( yn 1 )
implicit:
Some definitions
The method is said to be single step if it only uses the information from the last solution
point or within the current step to compute the new output value.
The Euler’s method or the Runga Kutta method are single step.
The method is said to be multiple step if uses information from several past solutions to
calculate the new output point.
The Adams method is a typical multiple step method.
Some definitions
Stiffness occurs in a problem where there are two or more very different scales of the
independent variable on which the dependent variables are changing.
In such cases even small value for the step size can cause instability.
dx
Exampl 1195x 1995 y, x(0) 2
e dt
dy
1197x 1997 y, y (0) 2
dt
The analytical solution x (t ) 10 e 2 t 8e 800t , y (t ) 6e 2 t 8e 800t
is
if we use Euler’s method to solve these xi 1 xi h f ( xi , y i ) xi 0.1(1195 xi 1995 yi )
yi 1 yi h g ( xi , yi ) yi 0.1(1197 xi 1997 yi )
This gives x(0.1)=640, y(0.1)=636
The exact values are x(0.1)=8.187, y(0.1)=4.912
Some definitions
dx dy
This equation group 1195x 1995 y 1197x 1997 y
dt dt
1195 1995
'
x x
A , where A
y y
can be 1197 1997
written
if we compute the eigenvalue for this matrix A they are 1= -2 and 2= -800.
When the eigenvalues of A have real parts that are negative and differ widely in magnitude
the system is stiff. As in this example.
The example above was linear. If we have a nonlinear x1 f1 ( x1 ...xn )
system
x2 f 2 ( x1 ...xn )
... ....
x
n f n ( x1 ...xn )
f i
we must consider the Jacobian matrix whose terms
x j
are
Example of an implicit method
for Stiff problems
What to do with our stiff problem? We can use an implicit
method.
The implicit form of the Euler’s method
xi 1 xi h f ( xi 1 , yi 1 ) xi 0.1(1195 xi 1 1995 yi 1 )
yi 1 yi h g ( xi 1 , yi 1 ) yi 0.1(1197 xi 1 1997 yi 1 )
We rewrite
1
xi 1 (1 1195(0.1) 1995(0.1) xi
y
i 1 1197 (1 1997(0.1) yi
We solve and get x(0.1)=8.23, y(0.1)=4.90
The exact values are x(0.1)=8.187,
y(0.1)=4.912.
We are reasonably close.
Runge-Kutta methods
Runge-Kutta methods are one of the most widely used methods.
The 4th order Runga-Kutta is the workhorse of ODE-solving.
In the fourth order Runga-Kutta the derivatives are evaluated 4 times at each step:
once at the initial point, twice at trial midpoints, and once at trial end point. From
these points the final function value is computed.
1
y n 1 yn (k1 2k 2 2k 3 k 4 )
6
k1 h f ( xn , y n )
1 1
k 2 h f ( xn h, yn k1 )
2 2
1 1
k 3 h f ( x n h, y n k 2 )
2 2
k 4 h f ( x n h, y n k 3 )
Example
dy
As an example: x y, y (0) 1, step size h 0.1
solve dx
k1 = 0.1*(0+1)=0.10000
k2 = 0.1*(0.05+1.05)=0.1100
k3 = 0.1*(0.05+1.055)=0.11050
k4 = 0.1*(0.10+1.11050)=0.12105
Y(0.1)=1.0+1/6*(0.1+0.22+0.221+0.12105)=
1.11034
Runge-Kutta
There are many variations of the Runga-Kutta
methods.
Runge-Kutta in a general form:
s
y n 1 y n bi k i
i 1
s
ki h f (tn ci h, yn aij k j )
j 1
by looking at parameters aij we can obtain 3 classes of RK
methods.
Explicit methods if aij = 0, j i
Semi-implicit methods if aij = 0, ji
Implicit methods No restrictions on aij
Runge-Kutta
Explicit methods are simple to use and do not require the solution of sets of non-linear
algebraic equations. They have a limited stability region and are thus suitable for non-
stiff problems.
Semi-implicit methods (SIRK) require iterative solution on stage calculations. Methods
possess very strong stability properties, which makes them useful for stiff problems. A
special case of these are the Diagonally Implicit RK methods (DIRK).
The implicit methods are attractive from the aspect of accuracy and stability but have
many implementation difficulties.
Multistep methods
Predictor-corrector methods store the solution along the way, and use those results to
extrapolate the solution one step advanced; they then correct the extrapolation using
derivative information at the new point.
According to the authors of Numerical Recipes: ”Only rather sophisticated predictor-
corrector methods are competitive”.
One class of linear multistep methods called ”the backward differentiation formulae”
BDF have been used succesfully in solving stiff problems. Implementations DIFSUB
(by C.W. Gear), GEAR (by A. Hindmarch) and FACSIMILE (by Curtis et al.).
Step size control
A ODE integrator should exert some adaptive control over its own progress, making
frequent changes in its stepsize.
Usually the purpose of the stepsize control is to achieve some predetermined accuracy in
the solution with minimum computational effort.
One idea for 4th order Runga-Kutta:
y ( x 2h ) y1 (2h )5 O (h 6 )
The errors for a step 2h and two times h
y ( x 2h ) y1 2(h )5 O (h 6 )
y2 y1 The s are scaled with h5, 0 is the desired
accuracy.
1/ 5
ho h1 0
1 This equation tells us how to increase or
decrease the stepsize for the next step.
How to Select the Method
Use Runge-Kutta when you do not know better.
When you need an easy implementation and you have full knowledge of you
problem and its properties you can use Euler.
For stiff problems use implicit methods (perhaps SIRK, BDF or DIRK).
If expensive right hand sides use multistep methods if your functions are
smooth.
Oscillatory, unstable, and discontinuous in t use single step methods.
Related docs
Other docs by hbf25307
Generating tangent linear and adjoint versions of NASAGMAO's Fortran - Download as PDF
Views: 4 | Downloads: 0
