Hjorth-Jensen M. Computational Physics

Подождите немного. Документ загружается.

16.2. DIFFUSION EQUATION 309

ENDDO

! write out for t = 0

t = 0 .

DO i =1 , ndim

WRITE(6 , ) t , i h , v ( i )

ENDDO

! setup the matrix to be inve rted

a = 0 . ; u =0.

DO i =1 , ndim 1

a ( i , i ) =1.+2. f a ctor

a ( i , i +1)= f a c t or

a ( i +1 , i )= f a c t or

ENDDO

a ( ndim , ndim ) = 1.+ 2. fa c t or

! now i n ver t the matrix

CALL matinv ( a , ndim , det )

DO i = 1 , m

DO l =1 , ndim

u( l ) = DOT_PRODUCT( a ( l , : ) ,v ( : ) )

ENDDO

v = u

t = i k

DO j =1 , ndim

WRITE(6 , ) t , j h , v ( j )

ENDDO

DEALLOCATE ( a ) ; DEALLOCATE ( u , v )

END SUBROUTINE solve_1dim_equation

END MODULE one_dim_heat_equation

PROGRAM heat_eq_1dm

USE one_dim_heat_equation

IMPLICIT NONE

INTERFACE

DOUBLE PRECISION FUNCTION f u n c t i o n _ i n i t i a l ( x )

IMPLICIT NONE

DOUBLE PRECISION, INTENT( IN) : : x

END FUNCTION f u n c t i o n _ i n i t i a l

END INTERFACE

CALL i n i t i a l i s e

OPEN(UNIT=6 , FILE=’ heat . dat ’ )

CALL solve_1dim_equation ( f u n c t i o n _ i n i t i a l )

310 CHAPTER 16. PARTIAL DIFFERENTIAL EQUATIONS

CLOSE(6)

END PROGRAM heat_eq_1dm

DOUBLE PRECISION FUNCTION f u n c t i o n _ i ni t i a l ( x )

IMPLICIT NONE

DOUBLE PRECISION, INTENT( IN) : : x

f u n c t i o n _ i n i t i a l = SIN ( x )

END FUNCTION f u n c t i o n _ i n i t i a l

16.2.4 Crank-Nicolson scheme

It is possible to combine the implicit and explicit methods in a slightly more general approach.

Introducing a parameter

(the so-called -rule) we can set up an equation

(16.46)

which for

yields the forward formula for the ﬁrst derivative and the explicit scheme, while

yields the backward formula and the implicit scheme. These two schemes are called the

backward and forward Euler schemes, respectively. For

we obtain a newscheme after its

inventors, Crank and Nicolson. This scheme yields a truncation in time which goes like

and it is stable for all possible combinations of and .

Using our previous deﬁnition of

we can rewrite the latter equation as

(16.47)

or in matrix-vector form as

(16.48)

where the vector

is the same as deﬁned in the implicit case while the matrix is

(16.49)

16.2.5 Non-linearterms and implementationofthe Crank-Nicoloson scheme

16.3 Laplace’s and Poisson’s equations

Laplace’s equation reads

(16.50)

16.3. LAPLACE’S AND POISSON’S EQUATIONS 311

with possible boundary conditions on the border. There is no time-dependence.

Choosing equally many steps in both directions we have a quadratic or rectangular grid, depend-

ing on whether we choose equal steps lengths or not in the

and the directions. Here we set

and obtain a discretized version

(16.51)

and

(16.52)

which we rewrite as

(16.53)

and

(16.54)

which gives when inserted in Laplace’s equation

(16.55)

This is our ﬁnal numerical scheme for solving Laplace’s equation. Poisson’s equation adds only

a minor complication to the above equation since in this case we have

and we need only to add a discretized version of resulting in

(16.56)

It is fairly straightforward to extend this equation to the three-dimensional case. Whether we

solve Eq. (16.55) or Eq. (16.56), the solution strategy remains the same. We know the values of

at or and at or but we cannot start at one of the boundaries

and work our way into and across the system since Eq. (16.55) requires the knowledge of at all

of the neighbouring points in order to calculate

at any given point.

The way we solve these equations is based on an iterative scheme called the relaxation

method. Its steps are rather simple. We start with an initial guess for

where all values

are known. To obtain a new solution we solve Eq. (16.55) or Eq. (16.56) in order to obtain a

new solution

. Most likely this solution will not be a solution to Eq. (16.55). This solution is

in turn used to obtain a new and improved

. We continue this process till we obtain a result

which satisﬁes some speciﬁc convergence criterion.

A simple example may help in visualizing this method. We consider a condensator with

parallel plates separated at a distance

resulting in e.g., the voltage differences

and . These are our boundary conditions and we

ask what is the voltage between the plates? To solve this problem numerically we provide

below a Fortran 90/95 program which solves iteratively Eq. (16.55).

312 CHAPTER 16. PARTIAL DIFFERENTIAL EQUATIONS

programs/chap16/program2.f90

! Program to solve the 2 dim Laplace equation using i t e r a t i o n .

! No time dependence .

! I n i t i a l c onditions are read in by the functi on i n i t i a l i s e

! such as number of st ep s in the x direction , y direction ,

! xmin and xmax , ymin and ymax . Here we employ a square l a t t i c e

! with equal number of steps in x and y dir e c tio n s

! Note the s t ru ctur e of th i s module , i t contains various

! subrou tines for i n i t i a l i s a t i o n of the problem and sol uti on

! of the PDE with a given i n i t i a l f unct io n for u ( x , y )

MODULE two_dim_laplace_equation

DOUBLE PRECISION, PRIVATE : : xmin , xmax , ymin , ymax

INTEGER, PRIVATE : : ndim , i t e r a t i o n s

DOUBLE PRECISION, ALLOCATABLE, DIMENSION( : , : ) , PRIVATE : : u ,

u_temp

CONTAINS

! t h is f unct ion reads in the s i ze of l a t t i c e , xmin , xmax ,

! ymin and ymax and the number of i t e r a t i o n s

SUBROUTINE i n i t i a l i s e

IMPLICIT NONE

WRITE( , ) ’ read in number of mesh points in x and y di re ction

’

READ( , ) ndim

WRITE( , ) ’ read in xmin and xmax ’

READ( , ) xmin , xmax

WRITE( , ) ’ read in ymin and ymax ’

READ( , ) ymin , ymax

WRITE( , ) ’ read in max number of i t e r a t i o n s ’

READ( , ) i t e r a t i o n s

END SUBROUTINE i n i t i a l i s e

SUBROUTINE solve_2dimlaplace_equation ( func )

DOUBLE PRECISION : : h , x , y , pi , length , d i f f

INTEGER : : i , j , l

INTERFACE

DOUBLE PRECISION FUNCTION func ( argument )

IMPLICIT NONE

DOUBLE PRECISION, INTENT( IN) : : argument

16.3. LAPLACE’S AND POISSON’S EQUATIONS 313

END FUNCTION func

END INTERFACE

! defi ne the step s i ze

h = ( xmax xmin ) /FLOAT( ndim +1)

length = xmax xmin

! al loc ate space for the vector u and the temporary vector to

! be upgraded in every i t e r a t i o n

ALLOCATE ( u ( ndim , ndim ) )

ALLOCATE ( u_temp ( ndim , ndim ) )

pi = ACOS( 1.)

! se t up of i n i t i a l conditions at t = 0 and boundary condi tions

u = 0 .

DO i =1 , ndim

x = i h pi / length

u( i , 1) = func ( x )

u( i , ndim ) = func ( x)

ENDDO

! i t e r a t i o n algorithm s t a r t s here

i t e r a t i o n s = 0

DO WHILE ( ( i t e r a t i o n s <= 20) .OR . ( d i f f > 0.00001) )

u_temp = u ; d i f f = 0 .

DO j = 2 , ndim 1

DO l = 2 , ndim 1

u( j , l ) = 0.2 5 ( u_temp ( j +1 , l )+u_temp ( j 1, l ) + &

u_temp ( j , l +1)+u_temp ( j , l 1) )

d i f f = d i f f + ABS( u_temp ( i , j ) u( i , j ) )

ENDDO

i t e r a t i o n s = i t e r a t i o n s + 1

d i f f = d i f f / ( ndim +1) 2

ENDDO

! write out r e s u l t s

DO j = 1 , ndim

DO l = 1 , ndim

WRITE(6 , ) j h , l h , u ( j , l )

ENDDO

DEALLOCATE ( u , u_temp )

END SUBROUTINE solve_2dimlaplace_equation

END MODULE two_dim_laplace_equation

PROGRAM laplace_eq_2dim

USE two_dim_laplace_equation

314 CHAPTER 16. PARTIAL DIFFERENTIAL EQUATIONS

IMPLICIT NONE

INTERFACE

DOUBLE PRECISION FUNCTION f u n c t i o n _ i n i t i a l ( x )

IMPLICIT NONE

DOUBLE PRECISION, INTENT( IN) : : x

END FUNCTION f u n c t i o n _ i n i t i a l

END INTERFACE

CALL i n i t i a l i s e

OPEN(UNIT=6 , FILE=’ lapla ce . dat ’)

CALL solve_2dimlaplace_equation ( f u n c t io n _ i n i t i a l )

CLOSE(6)

END PROGRAM laplace_eq_2dim

DOUBLE PRECISION FUNCTION f u n c t i o n _ i ni t i a l ( x )

IMPLICIT NONE

DOUBLE PRECISION, INTENT( IN) : : x

f u n c t i o n _ i n i t i a l = 100 SIN ( x )

END FUNCTION f u n c t i o n _ i n i t i a l

The important part of the algorithm is applied in the function which sets up the two-dimensional

Laplace equation. There we have a do-while statement which tests the difference between the

temporary vector and the solution

. Moreover, we have ﬁxed the number of iterations to be at

most . This is sufﬁcient for the above problem, but for more general applications you need to

test the convergence of the algorithm.

16.4 Wave equation in two dimensions

The

-dimensional wave equation reads

(16.57)

with

and we have assumed that we operate with dimensionless variables. Possible

boundary and initial conditions with are

(16.58)

16.4. WAVE EQUATION IN TWO DIMENSIONS 315

We discretize again time and position,

(16.59)

and

(16.60)

which we rewrite as

(16.61)

and

(16.62)

resulting in

(16.63)

If we assume that all values at times and are known, the only unknown variable

is and the last equation yields thus an explicit scheme for updating this quantity. Wehave

thus an explicit ﬁnite difference scheme for computing the wave function

. The only additional

complication in our case is the initial condition given by the ﬁrst derivative in time, namely

. The discretized version of this ﬁrst derivative is given by

(16.64)

and at

it reduces to

(16.65)

implying that . If we insert this condition in Eq. (16.63) we arrive at a special

formula for the ﬁrst time step

(16.66)

We need seemingly two different equations, one for the ﬁrst time step given by Eq. (16.66) and

one for all other time-steps given by Eq. (16.63). However, it sufﬁces to use Eq. (16.63) for all

times as long as we provide

using

(16.67)

in our setup of the initial conditions.

The situation is rather similar for the

-dimensional case, except that we now need to

discretize the spatial

-coordinate as well. Our equations will now depend on three variables

whose discretized versions are now

(16.68)

316 CHAPTER 16. PARTIAL DIFFERENTIAL EQUATIONS

and we will let and for the sake of simplicity. The equation with initial

and boundary conditions reads now

(16.69)

We have now the following discretized partial derivatives

(16.70)

and

(16.71)

and

(16.72)

which we merge into the discretized

-dimensional wave equation as

(16.73)

where again we have an explicit scheme with

as the only unknown quantity. It is easy to

account for different step lengths for and . The partial derivative is treated in much the same

way as for the one-dimensional case, except that we now have an additional index due to the

extra spatial dimension, viz., we need to compute

through

(16.74)

in our setup of the initial conditions.

16.4.1 Program for the

wave equation and applications

16.5 Inclusion of non-linear terms in the wave equation

Part II

Advanced topics

317