Download Coding and Numerical Integration of Ordinary Differential Equations (ODEs) and more Assignments Earth Sciences in PDF only on Docsity!
Chapter 11
Coding & Numerical Integration
of ODEs
11.1 Introduction
Translating algorithms into code takes art, patience, and hard work. The process can be greatly shortened by following a few rules to implement the algorithm. The sources of the algorithm are not of concern here although in practice we should be familiar with the intent and limitations of any algorithm, particularly when it relates to numerical computations. We will leave these issues for later as we are mainly concerned with programming.
11.2 Rules of thumbs
This section lists a few rules of thumbs to improve the coding of algorithms, and help in the homework problem.
11.2.1 Define the problem
Understand what is required in the computations, define the input, the output and the computations to be performed. Making sure the problem is well-defined goes a long way in determining what kind of input data and what sort of output is required.
11.2.2 Divide and conquer
Dividing the problem into functionally separate tasks helps in making large cod- ing projects manageable. Identify the repetitive tasks and consider isolating them into well-designed subroutines and functions. This subdivision is critical for the remainder of the programming tasks as it: divides into manageable units, allows
70 CHAPTER 11. CODING & NUMERICAL INTEGRATION OF ODES
validation to proceed in stages, and simplify code reuse. The details of the compu- tations will then be boxed into compartments and a higher level abstraction about the algorithm is then possible. Most numerical tasks follow the plan:
- data input
- perform computations
- output results
Each of these tasks may contain a long list of substasks but the divide and conquer strategy can be applied to those recursively.
11.2.3 Use top-down design
Start designing the problem from the top starting with the main program and filling in the subroutine as you go. This will help you clarify the flow of the code to yourself. Write shell routines in the initial go and start filling them and validating them in groups. Postpone the details of each routine as long as possible.
11.2.4 Use pseudo-code
Describe in pseudo code (i.e. programming language independent) the steps needed in the algorithms.
11.2.5 Decide on optimal data structure
In most numerical computations this will simply arrays as they are the most effi- cient data structures for many types of computations.
11.2.6 Conceptual simplicity
Avoid using obfuscated code even though it may seem great at the moment. Ob- fuscated code can be hard to read, and even the implementors would have a hard time remembering what it does 3 months after they coded it.
11.2.7 Validation
It is important to gain assurances that parts of the code are working properly in order to catch coding and (more ominously) algorithmic mistakes. It is then imperative that subroutines are validated as one is developing the code. If an error is introduced later it will useful to restrict its tracking to the untested portions of the code. Validation can take many forms:
72 CHAPTER 11. CODING & NUMERICAL INTEGRATION OF ODES
stages to predict the value of the unknown function at the next time level starting from the known initial conditions. All that is needed is the repeated evaluation of the function F~ on the right hand side of equation 11.1. The steps of the algorithms can be described as follows: Given the solution at w~n^ and the time step ∆t:
- Compute the gradient of the function at time tn and predict a first temporary value at tn + ∆ 2 t
F^ ~ (1)^ = F~ ( w~n, tn) , w~(1)^ = w~n^ + ∆t 2
F~ (1)^ (11.7)
- Compute the gradient of the function at time tn + ∆ 2 t using the previous temporary value of the solution w~(1)^ and compute a second temporary value at tn + ∆ 2 t
F^ ~ (2)^ = F~
( w ~(1), tn +
∆t 2
) , w~(2)^ = w~n^ +
∆t 2
F~ (2)^ (11.8)
- Compute the gradient of the function at time tn + ∆ 2 t using w~(2) and compute a first prediction of the solution at time tn+
F^ ~ (3)^ = F~
( w ~(2), tn +
∆t 2
) , w~(3)^ = w~n^ + ∆t F~ (3)^ (11.9)
- Compute the gradient of the function at time tn +∆t and compute the correct solution at time tn+1 using a linear combinations of the slopes computed so far
F^ ~ (4)^ = F~
( w ~(3), tn + ∆t
) , w~n+1^ = w~n^ + ∆t
F~ (1)^ + 2 F~ (2)^ + 2 F~ (3)^ + F~ (4)
11.5 Analysis
Here we will attempt to define the problem and analyzing the algorithm so as to start sketching out the design of our code.
11.6 Problem definition
We need to integrate the function w~ in time given its initial conditions w~ 0.
11.6. PROBLEM DEFINITION 73
11.6.1 Physical and Numerical Input Parameters
The parameters of the problem we wish to compute consist of
- the Coriolis parameter f
- the initial condition w~ 0
These parameters are intrinsic to the problem and do not depend on the numerical scheme. The numerical parameters consist of:
- The time step size ∆t
- and the number of time steps to take say M. The final integration time would then be tM = M ∗ ∆t assuming we start our timer at 0.
- We probably need to save the solution at intermediate steps to plot its vari- ation in time. We will thus need a parameter to determine the frequency of writing out the solution.
11.6.2 Output Parameters
The solution w~ need to be saved as the computations proceed. Diagnostics about the computations should be helpful in telling us whether we are proceeding on the right track or not. One possible diagnostics is the energy contained in the system. Multiplying the equations by u and v respectively and adding them up it is evident that
u
du dt
dv dt
d dt
( u^2 + v^2 2
) = 0 (11.11)
and hence the energy (u^2 + v^2 )/2 should remain constant at all times. It is thus a good idea to monitor the energy as the calculation proceeds.
11.6.3 Analysis of the RK4 algorithm
One can make the following remarks about the algorithm:
- The algorithm requires 4 steps that are similar in appearance.
- At each stage the slope is computed and a temporary value for the solution is computed.
- The final slope is a linear combinations of the slopes computed so far.
- The algorithm requires the repeated evaluation of the slope with different input parameters
11.7. THE CODING AND INCREMENTAL VALIDATION 75
routine Initialize routine save-solution routine diagnostics routine rk calculate rhs of equation and temporary value for each stage end routine rk routine rhs set physical parameter f find the slope for a given w and time routine rhs
11.7 The Coding and Incremental Validation
A first implementation of the main program, may look as follows
program inert implicit none integer, parameter :: nd=2! dimension of solution vector real :: w(nd)! solution vector integer :: Ntimes, isnap! number of time steps and save frequency
print *,’Enter dt, Ntimes, isnap’ read *, dt, Ntimes, isnap
call Initialize(w,nd) call PrintOut(w,nd) call Diagnostics(w,nd) do it = 1,Ntimes call Rk4(w,time,dt,nd) if (mod(it,isnap)==0) then call Printout(w,nd) call Diagnostics(w,nd) endif enddo
stop end program inert
There are a few variable declarations missing but the compiler will help us find those pretty easily because the implicit none will require that all variables be declared. Now we begin the task of writing the subroutines. A first cut at the initialize subroutine would look like:
76 CHAPTER 11. CODING & NUMERICAL INTEGRATION OF ODES
subroutine Initialize(w,nd) implicit none integer, intent(in) :: nd real, intent(out0 :: w(nd) w(1) = 1. w(2) = 0. return end subroutine Initialize
The subroutine is quite simple and it is unlikely that we have introduced a bug. Nevertheless, we could probably produce an executable and run the code to check it out. To do that however, we would need to provide the remainder of the sub- routines. At this point all we need is the declaration for the computer to produce the executables. These would look like:
!*********************************************************************** subroutine PrintOut(w,nd) implicit none integer, intent(in) :: nd real, intent(in) :: w(nd) return end subroutine PrintOut !*********************************************************************** subroutine Diagnostics(w,nd) implicit none integer, intent(in) :: nd real, intent(in) :: w(nd) return end subroutine Diagnostics !*********************************************************************** subroutine Rk4(w,time,dt,nd) implicit none integer, intent(in) :: nd! dimension of solution vector real, intent(in) :: time,dt! time level and time step real, intent(inout) :: w(nd)! solution vector return end subroutine Rk
These routines contain only variable declarations and no executable statetements at this stage. They are enough now to create an executable and run the code to initialize the vector w. At this point we can ascertain that the code runs but we have no way of verifying the output; we need to write out the solution and the routine PrintOut is then necessary. It is probably best if we write the solution to a file that we can examine after the run is done. Here it is:
78 CHAPTER 11. CODING & NUMERICAL INTEGRATION OF ODES
call rhs(r1,w,time,nd)! stage 1 wt = w + 0.5dtr call rhs(r2, wt, time+0.5dt,nd)! stage 2 wt = w + 0.5dtr call rhs(r3, wt, time+0.5dt,nd)! stage 3 wt = w + dtr call rhs(r4, wt, time+dt, nd)! stage 4 w = w + dt * ( r1 + 2.0(r2+r3) + r4 ) / 6.
return end subroutine Rk
Validating this code is relative simple provided we give it a trivial forcing func- tion. The solution should not change for example if the right hand side is zero. The right hand side subroutine for our case however would be as follows:
subroutine rhs(r,w,time,nd) implicit none integer, intent(in) :: nd real, intent(in) :: time, w(nd) real, intent(out) :: r(nd) real, parameter :: f = 1 r(1) = f * w(2) r(2) =-f * w(1) return end subroutine rhs
Notice that we chose to declare the Coriolis parameter locally since it is not needed by any of the other routines.
11.8 Run Plotting and Validation
Once the bugs and compiler errors are corrected, the code can be run for a short run. Initial runs should proceed with optimization parameters off, and debuggers parameters on. These should include
- -g to generate debug operation and inhibit optimization
- -C to check array bounds at run-time. The code stops with an informational message if the array bounds are trespassed.
- Turn floating point exception on to catch overflows, divide by zero, and other floating point excpetions. Underflow exception is not generally a serious problem and can be overlooked.
11.8. RUN PLOTTING AND VALIDATION 79
- Trap uninitialized variables. Some compilers have the options to insert a nonsensical values in all variables declared in the program. The code crashes if a variable is used before its value is set by the programmer. Unfortunately it does not seem like the Portland Group compiler has this feature while the Intel ”ifort” compiler does; It is -trapuv
- Another common source of errors is the mismatch between the dummy argu- ment list in a subroutine declaration and the calling statement. This can take the form of eiter wrong data type, wrong size, or wrong sequence. We will introduce an interface construct that help us catch these sort of problems at compile time.
Once the code has been validated for a small test run, the optimization can be removed since they tend to slow execution severely.
11.8.1 Plotting
The plotting in matlab is usually a breeze. Simply launch matlab and read in the 2 columns with the load matlab command. The two columns will contain the evolution of the solution caught at every isnap time steps.
