up
PDF

## Solving the advection PDE in explicit FTCS, Lax, Implicit FTCS and Crank-Nicolson methods for constant and varying speed. Accuracy, stability and software animation Report submitted for fulﬁllment of the Requirements for MAE 294 Masters degree project Supervisor: Dr Donald Dabdub, UCI.

### Contents

2.1.1 FTCS
2.1.2 Downwind
2.1.3 Upwind
2.1.4 LAX
2.1.5 Lax-Wendroﬀ
2.1.6 Leap-frog

2.2.1 Implicit FTCS
2.2.2 Wendrof
2.2.3 Crank-Nicolson

4.1 Case 1
4.2 Case 2
4.3 Case 3
4.4 Case 4
4.5 case 5
4.6 case 6
4.7 case 7
4.8 case 8

6.1 Plots
6.1.1 case 1
6.1.2 case 2
6.1.3 case 3
6.1.4 case 4
6.1.5 case 5
6.1.6 case 6
6.1.7 case 7
6.1.8 case 8
6.2 Source code

### 1 Introduction

The goal of this project is to analyze and compare diﬀerent numerical methods for solving the ﬁrst order advection PDE.

Following the analytical analysis for stability of the numerical scheme, animation were done to visually illustrate and conﬁrm these results. This was carried for diﬀerent parameters. The animation was programmed in Mathematica and saved to animated gif ﬁles which was then loaded into the HTML version of this report.

Fortran 95 was used for the computation part, while Mathematica was used for the animation and graphics part.

The above link contains all the supporting material for the project, including the Fortran program (in source and windows executable format) used to carry the main computation, and the Mathematica program used to do the animation and the Unix bash ﬁle used to process the computation for diﬀerent parameters.

The speciﬁc PDE example used for the analysis and animation was the one provided by Professor Donald Dabdub for the ﬁnal exam for his MAE 185 course (Numerical methods for mechanical engineers) in spring 2006. This PDE is described below:

Solve numerically

Where is the concentration of a given material as a function of time and space.

The above is solved for the following 2 cases

1.
(the advection speed, or the speed at which the mass is being transported) is a constant value given as (2 ft/min.
2.
is a function of time deﬁned as

The problem parameters are:

Where . Initial conditions are

The boundary conditions are

This PDE is an example of an IBVP (Initial and Boundary Value Problem).

Diﬀerent numerical methods are used to solve the above PDE. The methods are compared for stability using Von Neumann stability analysis.

The numerical methods are also compared for accuracy. This was done by comparing the numerical solution to the known analytical solution at each time step. The comparison was done by computing the root mean square error (RMSE) between the numerical and the analytical solution at each time step.

The method with the least RMSE at the end of the simulation is considered the most accurate.

The above PDE has a known analytical solution which is

The above analytical solution indicates that the initial concentration will move from left to right with the advection speed .

The formulation of each numerical method is shown below. is used to represent the space between 2 space grid point, or the space step size, and is used to represent the time step.

The space line has grid points. The spacing was ﬁxed at ft for all the methods and for all the test cases, while was changed. This made comparing the diﬀerent methods simpler. The following diagram illustrates the discretization used.

Should we consider the lower left and the lower right grid points above as part of the initial conditions, or part of the boundary conditions?

Stability of each method is derived. Stability is important, since by the Lax-Richtmyer equivalence theorem1 , stability implies convergence of the solution. Convergence of the numerical solutions implies that as the step size becomes smaller, the numerical solution converges to the analytical solution.

Explicit and implicit numerical methods are used. When solving for the future value of the solution at a single node in terms of only past values, the method is called an explicit method. In other words, when the only unknown is the future value of the solution at a single node, and everything else on the right hand side of the ﬁnite diﬀerence equation is a solution derived at earlier time step, the method is explicit.

An implicit method is one in which the ﬁnite diﬀerence equation contains the solution at a at future time at more than one node. In other words, future solution are being solved for at more than one node in terms of the solution at earlier time. Implicit methods therefor are usually solved by matrix methods by solving where represents present present known solution values, and are the unknown future solution values, and is the coeﬃcient matrix which will usually be block diagonal (or tri diagonal) in shape.

In the derivations below, the notation of is used to indicate the solution at time step and at space node Hence is written as . This notation seems to be more clear than the notation.

Diﬀerent ﬁnite diﬀerence schemes for solving a PDE are obtained by using diﬀerent methods of approximating the derivative terms in the PDE.

This will be illustrated using the space derivative This derivative can be approximated in one of the following 3 ways (all at time step )

#### 1.3 Center diﬀerence

The following are the derivation of a number of methods for solving the advection PDE obtained by using the above deﬁnitions for the derivative when applied to both space and time.

### 2 Numerical schemes

#### 2.1 Explicit Methods

##### 2.1.1 FTCS

With FTCS, the forward time derivative, and the centered space derivative are used. Hence the advection PDE can be written as

 (0)

Solving for results in

 (1)

This method will be shown to be unconditionally unstable.

##### 2.1.2 Downwind

Here, the forward time derivative for is used and also the forward space derivative for. This results in

This method will be shown to be unconditionally unstable as well.

##### 2.1.3 Upwind

Here, the forward time derivative for is used, and the backward derivative for is used. This results in

This will be shown to be stable if

##### 2.1.4 LAX

Looking at the FTCS eq (1) above, and shown below again

The term above is replaced by its average value to obtain the LAX method

 (4)

This method will be shown to be stable if

##### 2.1.5 Lax-Wendroﬀ

By using the second-order ﬁnite diﬀerence scheme for the time derivative, the method of Lax-Wendroﬀ method is obtained

##### 2.1.6 Leap-frog

In this method, the centered derivative is used for both time and space. This results in

This method requires a special starting procedure due to the term . Another scheme such as Lax can be used to kick start this method.

#### 2.2 Implicit Methods

##### 2.2.1 Implicit FTCS

Given the explicit FTCS derived above

The above is modiﬁed it by evaluating the space center derivative at time step instead of at time step , this results in

 (5A)

Hence

 (5B)

Writing it in matrix form, ﬁrst letting results in

Where is the number of space grid points.

The above is written as

Solving for which represents the solution at time step or at time . represents the current solution at time step , and is the matrix of the coeﬃcients shown above.

Due to the form of the matrix, (Called tri diagonal, or Block diagonal), an algorithm that takes advantages of this form is used. This is called the Thomas algorithm. This greatly speeds up the solution. If we had used a general algorithm to solve this system such as the Gauss elimination method, it would have been much slower, making the implicit method not practical to use. (Some tests on the same data showed the Thomas algorithm to be 50 times faster than Gaussian elimination).

##### 2.2.2 Wendrof

This method uses center diﬀerence for the derivative around the space step and the time step

This leads to the following scheme

This can also be solved using similar matrix method to that used for the implicit FTCS. This method is not used in this report.

##### 2.2.3 Crank-Nicolson

By taking the average of the explicit FTCS and the implicit FTCS formulations (shown again below), the C-N scheme is derived

Taking the average of the above results in

Now the system is setup to solve for future values as follows. Let , the system can be written as

Thomas algorithm is used to solve the above system for

### 3 Stability analysis

A numerical solution is stable if the ”energy content” remain below some limiting value no matter how long the solution is integrated. In essence, this means that the solution does not ’blow up’ after some time. This can be called BIBO stability (Bounded In Bounded Out).

Hence one way to analyze the stability of the numerical solution is to determine an expression that relates the amplitude of the solution between 2 time steps, and to determine if this ratio remain less than or equal to a unity as more and more time steps are taken.

This type of analysis is called Von Neumann stability analysis for numerical methods.

The analysis is based of ﬁnding an expression for the magniﬁcation factor of the wave amplitude at each step. The solution will be stable if this magniﬁcation factor is less than one.

Let the magniﬁcation factor be . The numerical scheme is stable iﬀ

The Courant–Friedrichs–Lewy (CFL) criteria for stability says that

Where and are as deﬁned above: is the wave speed, and

The number is also called the courant number .

Some numerical methods will be shown to be unconditionally unstable (such as explicit FTCS and the explicit upwind). This means that even if courant number was , the numerical solutions will eventually become unstable.

Some explicit methods such as LAX, are conditionally stable if the courant number was

Implicit methods are unconditionally stable, hence courant number is not used for these methods. However, this does not mean one can take as large step as one wants with the implicit methods, since accuracy will be aﬀected even if the solution remain stable.

Hence, the best numerical scheme is one in which the largest step size can be taken, with the least amount of inaccuracy in the numerical solution while remaining stable.

For numerical scheme that are conditionally stable, it can be seen from the CFL condition that for a ﬁxed speed and ﬁxed , the maximum time step that can be taken is given by

It can be immediately seen from above, that for the case when the advection speed is varying and is a function of time such as the case when implying that the speed is increasing with time, then when using a ﬁxed time step it will eventually become larger than and the numerical scheme will be unstable. This is because as is becoming larger and larger, while is ﬁxed, the term will become smaller and smaller.

Hence to keep the courant number , the time step taken must remain less than , hence using a ﬁxed time step with increasing will eventually lead to instability.

This will aﬀect the explicit methods that are conditionally stable such as the LAX method, since the Lax method is explicit and depends on satisfying the CFL all the time for its stability. Implicit methods are stable for any time step.

In the following we derive the details of the stability analysis and use Von Neumann analysis to derive an expression for the ampliﬁcation factor for diﬀerent numerical schemes.

So to summarize:

1.
Explicit FTCS is unconditionally unstable.
2.
Explicit LAX is stable if , or in other words,
3.
Implicit FTCS and C-R are stable for all

#### 3.1 Stability analysis for FTCS

Using Von Neumann method, the following trial solution to the PDE is assumed

where and is the wave number  and is the amplitude of the wave, as a function of time.

Hence the solution at time step and at is written as

 (2)

Substitute this trial solution (2) into the (1) results in

 (3)

Let be the ratio of the amplitude of the wave at time step relative to that at time step . hence

Divide (3) by results in

Divide the above by

Hence

This implies that regardless of the time step selected or the space step , hence

For a ﬁxed speed , the instability can be delayed by making smaller, but it could not be prevented. Eventually this numerical solution will blow up. This will be illustrated below in an animation. See case 3 and 4 as examples.

The instability can be delayed by making smaller for a ﬁxed or by making larger for a ﬁxed .

#### 3.2 Stability analysis of the downwind method

Substitute the trial solution into the above

Let

Hence

Hence for stability it is required that

or

since , a positive quantity, then the above condition can not be satisﬁed. Hence the downwind method is unconditionally unstable.

#### 3.3 Stability analysis of the upwind method

Substitute the trial solution into the above

Let

Hence

Hence

Hence for stability it is required that

or

Which will be true only if or hence this implies

Hence the upwind method is stable if the CFL condition is satisﬁed. This will be seen as the same stability condition for the Lax method below.

#### 3.4 Stability analysis of Lax

Replace the trial function from (2) in Lax formulation in (4) and obtain

Divide by , the magniﬁcation factor is obtained

Hence

Since and , then it is seen that if

Hence the following is the condition for stability

As mentioned earlier, this is called the CFL condition.

The Lax method is stable for however, a modiﬁed version of this method is more accurate, which is the Lax-Wendroﬀ method.

#### 3.5 Stability of Lax-Wendroﬀ

This is the same as the Lax method. The method is stable if

#### 3.6 Stability analysis of the Implicit FTCS

Replace the trial function from (2) in (5B) results in

Divide by

Hence

Hence this shows that the

This property is common to all implicit methods.

Even though the implicit FTCS is stable, it is not very accurate. See case 8 below for an example.

### 4 Solution Results and Output

For the Fortran implementation, the following methods are implemented. The explicit FTCS, Explicit Lax, Implicit FTCS, and Implicit Crank-Nicolson.

For each method, the following was generated

1.
CPU time used for the run.
2.
snap shot of the solution at and minutes.
3.
RMSE between the numerical solution and the analytical solution.
4.
Animation of the numerical solution. The animation was done by taking snapshots of the solution at regular intervals in Fortran. These were saved to disk. Then Mathematica was used to generate the animation and the plots.

To compare the stability and accuracy of the methods, the time step was changed (increased) and a new run was made. 8 diﬀerent values of time steps are used. So there are 8 tests cases. These 8 test cases were run for both ﬁxed speed and for ft/min.

This table below summarizes these cases.  The appendix contains all the plots. The animations are added as HTML links.

#### 4.1 Case 1

sec,

 Speed Method CPU time (sec) RMSE Animation (2D) plots U=2 Explicit FTCS 20 0.0546 HTML HTML Explicit LAX 31 0.0543 HTML Implicit FTCS 45 0.0548 HTML C-R 49 0.0544 HTML U=t/20 Explicit FTCS 21 0.003 HTML HTML Explicit LAX 31 0.0031 HTML Implicit FTCS 67 0.0031 HTML C-R 69 0.0032 HTML

Note the following: The explicit FTCS remained stable throughout the run due to the small time step. All other methods were stable as well during the run. For the CPU for the varying case, notice that for the implicit methods this value is larger than the CPU for the same methods but when is ﬁxed. This is due to the fact that the matrix is no longer constant, and must be recomputed at each time step before calling Thomas algorithm to solve system.

Also notice that the CPU time for the implicit methods is larger than the explicit methods. This is due to the extra computational cost in solving . Even when using Thomas algorithm, this is still more expensive than the explicit methods when number of time steps is large.

#### 4.2 Case 2

sec,

 Speed Method CPU time (sec) RMSE Animation (2D) plots U=2 Explicit FTCS 2.42 0.01264 HTML HTML Explicit LAX 3.48 0.0057 HTML Implicit FTCS 4.7 0.00742 HTML C-R 4.9 0.00575 HTML U=t/20 Explicit FTCS 2.5 0.00352 HTML HTML Explicit LAX 3.5 0.00329 HTML Implicit FTCS 7 0.00337 HTML C-R 7.5 0.0033 HTML

The explicit FTCS is stable for most of the run, near the end it is starting to be become unstable.

Notice that around 26 minutes that ”bubbles” are starting to show up in the numerical solution downstream. This is a characteristic of how this method becomes unstable.

This will be more clear in the next test cases when the time step is made larger. For the varying speed case, the explicit method using the same time step remained stable during the whole 30 minutes. This is because the average speed was less than 2 ft/min, hence the mass did not have to travel as long a distance as with ﬁxed speed of , and so the instability did not show up. Mathematically this can be explained by looking at the term , hence for smaller , the courant number is smaller. Notice also the RMSE is smaller for variable speed compared to ﬁxed speed. Again this is related to the smaller average speed making the courant number smaller.

#### 4.3 Case 3

In this case, we slightly make the time step longer than before. We start to see the instability of FTCS.

sec, ft for ﬁxed

 Speed Method CPU time (sec) RMSE Animation (2D) plots U=2 Explicit FTCS 1.9 0.0494 HTML HTML Explicit LAX 2.78 0.01125 HTML Implicit FTCS 3.7 0.01245 HTML C-R 3.9 0.01128 HTML U=t/20 Explicit FTCS 2.0 0.00365 HTML HTML Explicit LAX 2.9 0.00331 HTML Implicit FTCS 5.56 0.00346 HTML C-R 6 0.00331 HTML

For explicit FTCS, The solution now starting to show instability at 25 minutes. Lax remained stable since CFL is satisﬁed. Explicit FTCS is becoming less accurate as well. Explicit Lax is most accurate at this time step.

#### 4.4 Case 4

In this case, we slightly make the time step even longer than before. Now FTCS becomes more unstable.

sec, ft

 Speed Method CPU time (sec) RMSE Animation (2D) plots U=2 Explicit FTCS 1.73 0.15249 HTML HTML Explicit LAX 2.56 0.000563 HTML Implicit FTCS 3.34 0.009005 HTML C-R 3.45 0.00565 HTML U=t/20 Explicit FTCS 1.84 0.00380 HTML HTML Explicit LAX 2.53 0.00336 HTML Implicit FTCS 4,73 0.00358 HTML C-R 5 0.003373 HTML

FTCS Instability starts at around 20 minutes. LAX remained stable since CFL is satisﬁed. Lax remained the most accurate at this time step. It accuracy actually improved as the time step became larger.

#### 4.5 case 5

Again the time step is made longer than before. Now the explicit FTCS is completely unstable.

sec, ft

For the case of ﬁxed , we have , while for varying ,  the maximum value will be at the end of the run, which is ft/min., hence the CFL condition is changing, with a value of at the end of the run which is still

 Speed Method CPU time (sec) RMSE Animation (2D) plots U=2 Explicit FTCS 0.73 blows up HTML HTML Explicit LAX 0.281 0.000162 HTML Implicit FTCS 0.437 0.1306 HTML C-R 0.4 0.01028 HTML U=t/20 Explicit FTCS 0.28 blow up HTML HTML Explicit LAX 0.3 0.01117 HTML Implicit FTCS 0.40 0.0386 HTML C-R 0.4 0.01197 HTML

For the varying speed case, the explicit FTCS remained stable for the duration of the run as compared to the case with the ﬁxed speed. This is because the average wave speed is less than with the ﬁxed wave speed case.

The magniﬁcation factor depends on the speed of the wave.

With the varying speed case, the coeﬃcient was smaller during the whole run, since the maximum speed attained will be ft/min. as compared to ft/min. in the ﬁxed case.

We see than the smaller the speed the smaller the magniﬁcation (with everything else being ﬁxed).

If we have run the simulation a little longer for the varying speed case, we will see the instability with explicit FTCS. This below is a diagram showing 2 runs using the explicit FTCS both with ft/min, one was run for 30 minutes, and the second for 53 minutes. The run to 30 minutes showed no instability while the run for 53 minutes showed the instability. This show the explicit FTCS will eventually become unstable.

This is an animation of the above HTML

#### 4.6 case 6

In this case, the time step is increased so that is just above the CFL condition.

Notice now that the Explicit LAX method become unstable as expected. The other implicit methods remain stable. the explicit FTCS method now is completely unstable. The implicit FTCS method is starting to become less accurate.

sec, ft

 Speed Method CPU time (sec) RMSE Animation (2D) plots U=2 Explicit FTCS 0.7 blows up N/A blows up HTML Explicit LAX 0.25 0.1006 HTML Implicit FTCS 0.5 0.13945 HTML C-R 0.468 0.01104 HTML U=t/20 Explicit FTCS 0.28 blows up N/A blows up HTML Explicit LAX 0.31 0.04385 HTML Implicit FTCS 0.45 0.0428 HTML C-R 0.56 0.01317 HTML

Notice that explicit LAX takes much less CPU than any other method.

#### 4.7 case 7

sec, ft

 Speed Method CPU time (sec) RMSE Animation (2D) plots U=2 Explicit FTCS 0.65 blows up N/A blows up HTML Explicit LAX 0.9 blows up HTML Implicit FTCS 0.42 0.1531 HTML C-R 0.41 0.01244 HTML U=t/20 Explicit FTCS 0.265 blows up N/A blows up HTML Explicit LAX 0.29 0.01389 HTML Implicit FTCS 0.36 0.0493 HTML C-R 0.36 0.01525 HTML

Notice that the CPU for the implicit method when speed is ﬁxed is now higher than the CPU for the explicit methods. This can be explained as follows: since the time step now is larger than before, the number of times to solve has been reduced. This made the implicit methods faster.

This implies that

Using a relatively large time step, implicit methods become faster than the explicit methods.

#### 4.8 case 8

sec, ft

 Speed Method CPU time (sec) RMSE Animation (2D) plots U=2 Explicit FTCS 0.5 blows up N/A blows up HTML Explicit LAX 0.89 blows up HTML Implicit FTCS 0.453 0.1653 HTML C-R 0.36 0.01403 HTML U=t/20 Explicit FTCS 0.234 blows up N/A blows up HTML Explicit LAX 0.2187 0.01564 HTML Implicit FTCS 0.344 0.0557 HTML C-R 0.312 0.0174 HTML

#### 4.9 CPU comparison tables

As expected, CPU time usage will be less as the time step is increased. There is an anomaly cased noticed where the CPU time increased for the Lax method when the time step is increased from to This needs further investigation.

This table below summarizes the CPU time in seconds used by each method for the case of constant speed as time step is increased.

This is the plot of the above table

This table below summarizes the CPU time in seconds used by each method for the case of varying speed as time step is increased.

This is the plot of the above table

This plot below compares the CPU time for each method when the speed is constant vs. when the speed was changing with time.

#### 4.10 Accuracy comparison tables

This table below summarizes the RMS error from each numerical method as a function of changing the time step size. This is for case of constant speed.

Notice that the Lax method became more accurate when the time step was increased from to seconds, then it starts to become less accurate as time step is increased. This is counter intuitive to what one can expect. It will be interesting to investigate this further to obtain a mathematical explanation for this strange phenomena.

The accuracy of the implicit FTCS, and C-R also increased slightly as the time step became larger from to , then the implicit FTCS became worst in terms of accuracy as the time step increased.

C-R method accuracy did not deteriorate as much with increasing the time step. This shows the C-R scheme to be more robust.

This table below summarizes the RMS error from each numerical method as a function of changing the time step size. This is for case of changing speed.

The eﬀect of having the speed deﬁned as is to delay instability for the explicit methods as time step is increased. Notice also here the case where the Lax method became more accurate as the time step is increased from to .

### 5 Conclusion

4 diﬀerent numerical ﬁnite diﬀerence schemes are examined for CPU time, stability and accuracy in solving the advection PDE for constant speed and for a speed which is a function of time.

For accuracy, an interesting result is observed. The Lax scheme is the most accurate for Courant number close to unity. This means as the time step is increased, the Lax become more accurate of the 4 methods. But beyond the CFL condition, Both explicit methods (FTCS and Lax) became less accurate. Explicit FTCS became unstable sooner than Lax, while the implicit methods remained stable.

The implicit FTCS was less accurate than the C-R method.  This implies that one should use the Lax method if one can be satisﬁed with a time step such that the courant number is close to a unit.

For stability, Crank-Nicolson was the most stable of all methods. Stability by itself is not suﬃcient condition to use to select a numerical scheme. It must also be accurate. The C-R method has both these properties for the range of the time steps considered. But as mentioned above, there is a range of time steps in which the Lax method is more accurate than all the other methods.

For CPU usage, the explicit methods used less CPU time when the time step was small, up to This can be explained as follows: for small step size,  the number of time to solve  is large. Hence the implicit methods will be slower. As the time step is increased to the range of and over, the implicit methods actually became more CPU eﬃcient due to the fact that the number of times to solve is less because the number of steps is less.

In conclusion, the selection of a ﬁnite diﬀerence scheme depends on many factors. Stability and accuracy being the most important. The time step size plays a critical rule. For Courant number close to a unity, the Lax method is the most attractive. For larger time steps, the C-R method should be considered.

### 6 Appendix

#### 6.2 Source code

!*******************************************
!*
!* Solve the advection PDE using Explicit FTCS,
!* Explicit Lax, Implicit FTCS, and implicit Crank-Nicolson
!* methods for constant and varying speed.
!*
!*  Solve  dc/dt = -u dc/dx
!*  u = t/20 ft/minute
!*  and
!*  u constant
!*
!* Compiler used: gnu 95 (g95) on Cygwin. Gcc 3.4.4
!* Date: June 20 2006
!*
!* by Nasser Abbasi
!*******************************************

IMPLICIT NONE

REAL     :: DT,DX,max_run_time,length,snapshot_delta,
first_limit,second_limit
INTEGER  :: N,SNAPSHOTS
character(10) :: cmd_arg ! to read time step from command line

INTEGER :: method  ! 1=FTCS, 2=LAX, 3=Implicit FTCS, 4=C-R
INTEGER :: mode    ! 1=Fixed wind speed, 2=speed function of time

REAL    :: t_start, t_end, cpu_time_used,end_line(1002)
INTEGER :: ALL_DATA_FILE_ID
PARAMETER(ALL_DATA_FILE_ID=900)

! Initialize data. All methods will use the same
! parameters to make comparing them easier

! read delta t from command line.
CALL getarg(1,cmd_arg)
cmd_arg=TRIM(cmd_arg)
print *,= , cmd_arg
read(cmd_arg,*)dt   !delta in time,  in minutes

print *,Dt=,DT

N      = 1000      ! number of grid points in space
length = 100       ! length of space solution in feet

first_limit  = 0.25*length
second_limit = 0.35*length

DX = length/N     ! delta in space, in feets

max_run_time = 30.0  ! how long to run for in minutes
SNAPSHOTS    = 200   ! number of snapshots per run. Used for animation

snapshot_delta = max_run_time / SNAPSHOTS   ! time between each snap shot

print *,DT=,DT, minutes, DX=,DX, feets
print *,taking snapshots every , snapshot_delta , minutes

DO mode=1,2
print*,=======> processing mode ,mode
DO method=1,4   ! No enumeration data types in Fotran 90

CALL CPU_TIME(t_start)  ! get current CPU time
CALL process(mode,method,N,DT,DX,max_run_time,snapshot_delta,
first_limit,second_limit)
CALL CPU_TIME(t_end)    ! get current CPU time

cpu_time_used = t_end - t_start

WRITE(*,FMT=(A,I2,A,F12.5)) CPU TIME used for method, method,  = , cpu_time_used
! Now record test case parameters in last line
end_line=0
end_line(1)=cpu_time_used
end_line(2)=DT
end_line(3)=DX
end_line(4)=mode
end_line(5)=method

WRITE(UNIT=ALL_DATA_FILE_ID,FMT=*) end_line
CLOSE(ALL_DATA_FILE_ID)

END DO
END DO

!************************************
!*
!*
!************************************
SUBROUTINE process(mode,method,N,DT,DX,max_run_time,snapshot_delta,
first_limit,second_limit)
IMPLICIT NONE

INTEGER, INTENT(IN) :: mode,method,N
REAL,    INTENT(IN) :: DT,DX,max_run_time,snapshot_delta,
first_limit,second_limit

INTEGER :: I
LOGICAL :: snap_shot_at_15_taken
INTEGER :: ALL_DATA_FILE_ID
PARAMETER(ALL_DATA_FILE_ID=900)
REAL    :: snap_current_time
REAL    :: current_time
REAL    :: C(N)       ! current solution
REAL    :: CNEW(N)    ! future solution
REAL    :: CEXACT(N)  ! current exact solution
REAL    :: current_first_limit
REAL    :: A(N,N),aa(N),b(2:N),cc(N-1),CTEMP(N)  ! for C-R and implicit FTCS
REAL    :: K,speed
REAL    :: error,RMS   ! root mean square error between current and initial sol.

current_time      = 0.
snap_current_time = 0.

CALL initialize_solution(C,N,DX,first_limit,second_limit)
CEXACT = C
current_first_limit = first_limit

CALL pre_loop_initialization(mode,method,current_time,K,
DT,DX,N,C,ALL_DATA_FILE_ID,
A,aa,b,cc )

snap_shot_at_15_taken=.FALSE.

DO WHILE(current_time < max_run_time)

IF( snap_current_time >= snapshot_delta ) THEN
snap_current_time = 0.
WRITE( UNIT=ALL_DATA_FILE_ID, FMT=*) current_time, error, C
END IF

SELECT CASE(method)

CASE( 1:2 )

IF(method==1) THEN  ! ftcs
IF(mode==2)THEN
K  = speed(mode,current_time)*DT/(2.*DX)
ENDIF

DO I = 2,N-1
CNEW(I) = C(I) - K * ( C(I+1) - C(I-1) )
END DO
ELSE    !lax
IF(mode == 2) THEN
K = speed(mode,current_time)*DT/(DX)
ENDIF

DO I = 2,N-1
CNEW(I) = C(I) - K/2. * ( C(I+1) - C(I-1) ) +
(K**2.)/2 * ( C(I+1) +C(I-1)-2.*C(I) )
END DO
END IF

CNEW(1) = C(1)
CNEW(N) = C(N)  ! Boundary conditions
C=CNEW

CASE( 3 )  ! implicit ftcs

IF( mode == 2) THEN ! only need to update Matrix for varying U
K  = speed(mode,current_time)*DT/(2.*DX)

CALL init_A_matrix(A,K,N)
CALL init_diagonal_vectors(N,A,cc,aa,b)
END IF

CALL solve_thomas_algorithm(N,aa,b,cc,C,CNEW)
C = CNEW

CASE( 4 )   ! C-R

IF(mode == 2) THEN  !only need to update A if U changes
K  = speed(mode,current_time)*DT/(4*DX)   ! C-R
CALL init_A_matrix(A,K,N)
CALL init_diagonal_vectors(N,A,cc,aa,b)
END IF

CTEMP(1) = C(1)
CTEMP(N) = C(N)

DO I=2,N-1
CTEMP(I)=C(I)+K*C(I-1)-K*C(I+1)
END DO

CALL solve_thomas_algorithm(N,aa,b,cc,CTEMP,C)

END SELECT

IF( current_time>=15.0 .AND. (.NOT. snap_shot_at_15_taken)) THEN
snap_shot_at_15_taken = .TRUE.
CALL take_one_snap_shot(mode,method,15,N,C,DX)
END IF

current_time = current_time + DT
current_first_limit = current_first_limit + speed(mode,current_time)*DT
CALL get_current_exact_solution(CEXACT,N,current_first_limit,DX)
error = RMS(CEXACT,C,N)

snap_current_time = snap_current_time + DT

END DO

CALL take_one_snap_shot(mode,method,30,N,C,DX)

END SUBROUTINE process
!************************************
!*
!*
!************************************
SUBROUTINE pre_loop_initialization(mode,method,current_time,K,
DT,DX,N,C,ALL_DATA_FILE_ID,
A,aa,b,cc)
IMPLICIT NONE

INTEGER, INTENT(IN)  :: mode,method,N,ALL_DATA_FILE_ID
REAL,    INTENT(IN)  :: C(N),DT,DX,current_time
REAL,    INTENT(OUT) :: K,A(N,N),aa(N),b(2:N),cc(N-1)
REAL                 :: speed

SELECT CASE(method)
CASE( 1 )   ! FTCS

K  = speed(mode,current_time)*DT/(2.*DX)

IF(mode==1) THEN
OPEN(UNIT=ALL_DATA_FILE_ID, file=expAll.txt) ! all time shots
CALL print_to_file(C,exp0.txt,N,DX)
ELSE
OPEN(UNIT=ALL_DATA_FILE_ID, file=exp_extraAll.txt) ! all time shots
CALL print_to_file(C,exp_extra0.txt,N,DX)
END IF

CASE( 2 )    ! Lax

K = speed(mode,current_time)*DT/(DX)

IF(mode==1) THEN
OPEN(UNIT=ALL_DATA_FILE_ID, file=laxAll.txt) ! all time shots
CALL print_to_file(C,lax0.txt,N,DX)
ELSE
OPEN(UNIT=ALL_DATA_FILE_ID, file=lax_extraAll.txt) ! all time shots
CALL print_to_file(C,lax_extra0.txt,N,DX)
END IF

CASE( 3 )    ! Implicit FTCS

K  = speed(mode,current_time)*DT/(2.*DX)

CALL init_A_matrix(A,K,N)
CALL init_diagonal_vectors(N,A,cc,aa,b)

IF(mode==1) THEN
OPEN(UNIT=ALL_DATA_FILE_ID, file=impAll.txt) ! all time shots
CALL print_to_file(C,imp0.txt,N,DX)
ELSE
OPEN(UNIT=ALL_DATA_FILE_ID, file=imp_extraAll.txt) ! all time shots
CALL print_to_file(C,imp_extra0.txt,N,DX)
END IF

CASE( 4 )     ! C-R

K  = speed(mode,current_time)*DT/(4*DX)   ! C-R

CALL init_A_matrix(A,K,N)
CALL init_diagonal_vectors(N,A,cc,aa,b)

IF(mode==1) THEN
OPEN(UNIT=ALL_DATA_FILE_ID, file=crAll.txt) ! all time shots
CALL print_to_file(C,cr0.txt,N,DX)
ELSE
OPEN(UNIT=ALL_DATA_FILE_ID, file=cr_extraAll.txt) ! all time shots
CALL print_to_file(C,cr_extra0.txt,N,DX)
END IF
END SELECT

WRITE( UNIT=ALL_DATA_FILE_ID, FMT=*) current_time,0, C

END SUBROUTINE pre_loop_initialization
!************************************
!*
!*
!************************************
SUBROUTINE init_diagonal_vectors(N,A,cc,aa,b)
IMPLICIT NONE

INTEGER, INTENT(IN) ::N
REAL, INTENT(IN)    ::A(N,N)
REAL, INTENT(OUT)   ::aa(N),b(2:N),cc(N-1)

INTEGER ::I,J

J=2
DO I=1,N-1
cc(I)=A(I,J)
J=J+1
END DO
cc(1)=0

DO I=1,N
aa(I)=A(I,I)
END DO

J=1
DO I=2,N
b(I)=A(I,J)
J=J+1
END DO

END SUBROUTINE init_diagonal_vectors
!************************************
!*
!*
!************************************
SUBROUTINE initialize_solution(C,N,DX,first_limit,second_limit)
IMPLICIT NONE

INTEGER, INTENT(IN)    :: N
REAL,    INTENT(IN)    :: DX,first_limit,second_limit
REAL,    INTENT(INOUT) :: C(0:N-1)

INTEGER :: I
REAL    :: x, PI,av,R

PARAMETER( PI = ACOS(-1.) )

x = 0
av = (second_limit+first_limit)/2.0
R  = av - first_limit

C = 0.0

DO I=0,N-1

IF( x >= first_limit .AND. x <= second_limit ) THEN
C(I) = 1 + COS( PI * (x-av)/R )
END IF

x = x + DX
END DO

END SUBROUTINE initialize_solution
!************************************
!*
!*
!************************************
SUBROUTINE print_to_file(C,file_name,N,DX)
IMPLICIT NONE

REAL,    INTENT(IN) :: C(N),DX
INTEGER, INTENT(IN) :: N

CHARACTER* (*), INTENT(IN) :: file_name

INTEGER :: I
INTEGER :: FILE_ID
PARAMETER(FILE_ID=999)
REAL :: current_position

OPEN(UNIT=FILE_ID, file=file_name)

current_position = 0;
DO I=1,N

WRITE( UNIT=FILE_ID, FMT=* ) current_position ,\t, C(I)
current_position = current_position + DX

END DO

CLOSE(FILE_ID)

END SUBROUTINE print_to_file
!************************************
!*
!*
!************************************
SUBROUTINE init_A_matrix(A,K,N)
IMPLICIT NONE

INTEGER, INTENT(IN) ::N
REAL, INTENT(IN)    ::K
REAL, INTENT(OUT)   ::A(N,N)

INTEGER ::I

DO I = 2,N-1
A(I,I-1)   = -K
A(I,I) = 1
A(I,I+1) = K
END DO

A(1,1) = 1
A(N,N) = 1

END SUBROUTINE init_A_matrix
!************************************
!*
!*
!************************************
SUBROUTINE solve_thomas_algorithm(N,aa,b,c,old_c,new_c)
IMPLICIT NONE

REAL, INTENT(IN)    :: aa(N),b(2:N),c(N-1),old_c(N)
INTEGER, INTENT(IN) :: N
REAL, INTENT(INOUT) :: new_c(N)

INTEGER  :: I
REAL :: alpha(N),beta(2:N),g(N)

alpha(1) = aa(1)
DO I=2,N
beta(I)=b(I)/alpha(I-1)
alpha(I)=aa(I)-beta(I)*c(I-1)
END DO

g(1)=old_c(1)
DO I=2,N
g(I)=old_c(I)-beta(I)*g(I-1)
END DO

new_c(N)=g(N)/alpha(N)
DO I=N-1,1,-1
new_c(I)=(g(I)-c(I)*new_c(I+1))/alpha(I)
END DO

END SUBROUTINE solve_thomas_algorithm
!************************************
!*
!*
!************************************
REAL FUNCTION speed(MODE,time)
IMPLICIT NONE

INTEGER, INTENT(IN) :: MODE
REAL,    INTENT(IN) :: time

IF( MODE == 1 ) THEN
speed=2.0
ELSE
speed=time/20.0
END IF

END FUNCTION speed
!************************************
!*
!*
!************************************
SUBROUTINE take_one_snap_shot(mode,method,TIME,N,C,DX)
IMPLICIT NONE

INTEGER, INTENT(IN) ::TIME,mode,method,N
REAL,    INTENT(IN) ::C(N),DX

IF(TIME==15) THEN
SELECT CASE(method)
CASE(1)
IF(mode==1) THEN
CALL print_to_file(C,exp15.txt,N,DX)
ELSE
CALL print_to_file(C,exp_extra15.txt,N,DX)
END IF
CASE(2)
IF(mode==1) THEN
CALL print_to_file(C,lax15.txt,N,DX)
ELSE
CALL print_to_file(C,lax_extra15.txt,N,DX)
ENDIF
CASE(3)
IF(mode==1) THEN
CALL print_to_file(C,imp15.txt,N,DX)
ELSE
CALL print_to_file(C,imp_extra15.txt,N,DX)
END IF
CASE(4)
IF(mode==1) THEN
CALL print_to_file(C,cr15.txt,N,DX)
ELSE
CALL print_to_file(C,cr_extra15.txt,N,DX)
END IF
END SELECT
ELSE
SELECT CASE(method)
CASE(1)
IF(mode==1) THEN
CALL print_to_file(C,exp30.txt,N,DX)
ELSE
CALL print_to_file(C,exp_extra30.txt,N,DX)
END IF
CASE(2)
IF(mode==1) THEN
CALL print_to_file(C,lax30.txt,N,DX)
ELSE
CALL print_to_file(C,lax_extra30.txt,N,DX)
ENDIF
CASE(3)
IF(mode==1) THEN
CALL print_to_file(C,imp30.txt,N,DX)
ELSE
CALL print_to_file(C,imp_extra30.txt,N,DX)
END IF
CASE(4)
IF(mode==1) THEN
CALL print_to_file(C,cr30.txt,N,DX)
ELSE
CALL print_to_file(C,cr_extra30.txt,N,DX)
END IF
END SELECT
END IF
END SUBROUTINE    take_one_snap_shot
!************************************
!*
!*
!************************************
REAL FUNCTION RMS(CEXACT,C,N)
IMPLICIT NONE

REAL, INTENT(IN)    :: CEXACT(N),C(N)
INTEGER, INTENT(IN) :: N

INTEGER  :: I

RMS=0.
DO I=1,N
RMS = RMS+(CEXACT(I)-C(I))**2
END DO

RMS = RMS/N
RMS = SQRT(RMS)
END FUNCTION RMS
!************************************
!*
!*
!************************************
SUBROUTINE get_current_exact_solution(CEXACT,N,current_first_limit,DX)
IMPLICIT NONE
REAL, INTENT(IN)    :: current_first_limit,DX
REAL, INTENT(OUT)   :: CEXACT(0:N-1)
INTEGER, INTENT(IN) :: N

INTEGER :: I
REAL :: first_limit
REAL :: second_limit
REAL :: av,R,shift,x,PI

PARAMETER( PI = ACOS(-1.) )

first_limit  = 25.0
second_limit = 35.0

shift = current_first_limit - FIRST_LIMIT
first_limit = current_first_limit
second_limit = second_limit + shift

av = (second_limit+first_limit)/2.0
R  = av - first_limit

CEXACT = 0.
x = 0.
DO I = 0,N-1

IF( x >= first_limit .AND. x <= second_limit ) THEN
CEXACT(I) = 1 + COS( PI * (x -av)/R )
END IF

x = x + DX
END DO
END SUBROUTINE get_current_exact_solution

### 7 References

1.
Numerical Methods for physics. Second edition. Alejandro Garcia
2.
Applied Numerical Methods for Engineers. Terrence Akal.
3.
Computational Techniques for ﬂuid dynamics. Second edition. C.A.J.Fletcher