浏览： 日期：2020-06-10

M 3/4/5 SC – SCIENTIFIC COMPUTING -- Due: 3 May 2013

Exercise #3 1 11 March 2013

Exercise #3 Part – II due 22:00pm Friday 3 May

(worth 15 % for M3SC and 40% for M 4/5 SC)

Background and discussion

The discrete Real Fourier series represents the periodic real data xj at N discrete points as the sum

of N cosine and sine functions in the form :

(II-1)

This would be a N/2-size discrete Complex Fourier Series if yk = ak + I bk . 0 < k < N/2

Although there are purely real algorithms to calculate xj given ak and bk and vice versa, most

commonly these problems are mapped into a procedure utilizing the Discrete Complex Fourier

Series and Transform.

One of these mappings take two real periodic data sets x1j and x2j 0 < j < N and packs them

into a periodic complex data set zj = x1j + x2j · I of size N and performs the size N Discrete

Complex Fourier Series on this data set: t = CN z . The resultant complex t vector can be

unpacked into the two real size N/2 Fourier Series coefficients by exploiting the fact that for a real

data series xj , the resultant complex yk have the property: yk = y*N-k . [II-2]

For a purely imaginary set of data, the resultant complex yk have the property: yk = - y*N-k .

Thus the yk for each real data set can be recovered for 0 < k < N/2 . In detail:

y 1k = ( tk + t*N-k )/N and y 2k = -I · ( tk - t*N-k )/N 0 < k < N/2.

To recover x1j and x2j 0 < j < N , given y1k and y2k 0 < k < N/2, extend y1k and y2k from

0 < k < N/2 to 0 < k < N using relation [II-2] and then create the complex vector

wk = y1k + y2k · I for 0 < k < N and perform the Discrete Fourier Transform CN-1 w .

The real and imaginary parts of the resultant complex array are the two desired real periodic data.

8. Write the C function to calculate the Discrete Fourier Series of two real data sets, each of

size N and find the two complex vectors z1 and z2 of size N/2 of their Fourier coefficients.

This function should have the prototype:

int TwoReal2Complex(double *R1, double *R2,

complex double *Z1,complex double *Z2,

complex double *w, complex double *Wp, int N)

where: *R1 is a pointer to the type double array to hold R1 from R1[0] to R1[N-1];

*R2 is a pointer to the type double array holding R2 from R2[0] to R2[N-1];

*Z1 is a pointer to the type complex double array holding Z1 from Z1[0] to Z1[N/2];

*Z2 is a pointer to the type complex double array holding Z2 from Z2[0] to Z2[N/2];

*w is a pointer to the type double array w[0] to w[your choice!] that can be used as

vector temporary storage for intermediate complex works required by your

FastCS(… function.

*Wp is a pointer to the vector or powers of WN created in the section 1 function above

N is the size of the desired transform.)

M 3/4/5 SC – SCIENTIFIC COMPUTING -- Due: 3 May 2013

Exercise #3 2 11 March 2013

9. Write the C function to calculate the Discrete Fourier Transform of two real periodic data sets,

each of size N given the two complex vectors z1 and z2 of size N/2 of their Fourier

coefficients This function should have the prototype:

int Complex2TwoReal(complex double *Z1,complex double *Z2,

double *R1, double *R2, complex double *w,

complex double *Wp, int N)

where

*Z1 is a pointer to the type complex double array holding Z1 from Z1[0] to Z1[N/2];

*Z2 is a pointer to the type complex double array holding Z2 from Z2[0] to Z2[N/2];

*R1 is a pointer to the type double array to hold R1 from R1[0] to R1[N-1];

*R2 is a pointer to the type double array holding R2 from R2[0] to R2[N-1];

*w is a pointer to the type double array w[0] to w[your choice!] that can be used as

vector temporary storage for intermediate complex works required by your

FastCS(… function.

*Wp is a pointer to the vector or powers of WN created in the section 1 function above

N is the size of the desired transform.)

10. Poisson’s Equation: 2 = - describes many fields in physics: electric field potential just

one example. A distribution of charge generates an electric field in which the charges

move in obeyance to Gauss’s laws (actually, Maxwell’s Laws) .

In Cartesian (x,y,z) co-ordinates this equation becomes

2

2

2

2

2

2

x

y

z

. (10.1)

If we write the vector r = x i + y j + z k and the vector K = kx i + ky j + kz k where i, j and k

are the usual unit vectors we can solve equation (8) by the separation of variables technique.

Let us assume that (x,y,z) = ( r) = exp( I K · r ) ( K) for some K vector. And assume that

(x,y,z) = ( r) = exp( I K · r ) ( K) . Substituting these expressions into

equation (8) leads to the relation

- (kx · kx + ky · ky + kz · kz) ( K) = - ( K) (10.2)

for every K vector.

Assume we want to solve equation (8.1) in a triply-periodic unit cube. This restricts kx, ky and kz

to take the discrete values: 2 π m, 2 π n and 2 π p for the integers m, n & p.

(We omit using o as an integer label as it may become confused with zero!)

Let (x,y,z) and (x,y,z) be defined at the vertices of a periodic cubical grid of size 1 with an

internal spacing of size = (1/N). So our functions and are defined at the points where

(x,y,z) take the values ( l/N, m/N, n/N), We can label each vertex value lmn and lmn.

This allows us to write : as a sum of periodic modes in each space dimension:

(x,y,z) = ( j, k , l ) =

1

, ,

0

N

mnp

m n p

exp(? ? ? ? ?

?

) exp(? ? ? ? ?

?

) exp(? ? ? ? ?

?

)

and (10.3)

(x,y,z) = ( j , k , l ) =

1

, ,

0

N

mnp

m n p

exp(? ? ? ? ?

?

) exp(? ? ? ? ?

?

) exp(? ? ? ? ?

?

)

for 0 < j, k, l < N .We assume a given periodic charge is assigned a value jkl at each vertex.

The problem is to solve Poisson’s equation (8.1) to find jkl at each (x,y,z) grid point. Both and

are assumed to have a zero mean value in our unit cube and to be real everywhere.

Find the electric potential ( l/N, m/N, n/N) for the charge distribution:

= -200 3/8 < x < 1/2 , 3/8 < y < 5/8, 7/16 < z < 9/16,

M 3/4/5 SC – SCIENTIFIC COMPUTING -- Due: 3 May 2013

Exercise #3 3 11 March 2013

= 0 x=0

= 200 1/2 < x < 5/8 , 3/8 < y < 5/8, 7/16 < z < 9/16,

You may assume that =0 everywhere else in the cube.

Find the solutions jkl for N = 16, 32, 48, 64, 80, 96,128, 160, 192, 256, 320 … until the problem

is too big for the PC.

Suggestions:

a. Allocate type double arrays of size N3 to hold the input, the output and the work arrays.

{I would recommend allocating these three arrays as single dimension arrays of size N3 .

Use the mapping p = j+ Nk + N2l to find the location of vertex j,k,l }

b. Set the jkl array to zero. For j from 3N/8 to N/2-1, for k from 3N/8 to 5N/8 and for

l from 7N/16 to 9N/16 set jkl = -200. Repeat for j from N/2+1 to 5N/8 with jkl = 200.

c. Find the Discrete Fourier transform of this distribution jkl by:

i. For k from 0 to N-1, For l from 0 to N-1, do the size N Discrete Fourier Transform

starting at 0kl with an initial skip of 1.

ii. For j from 0 to N-1, For l from 0 to N-1, do the size N Discrete Fourier Transform

starting at j0l with an initial skip of N.

iii. For j from 0 to N, For k from 0 to N, do the size N Discrete Fourier Transform

starting at jk0 with an initial skip of N2.

d. Calculate mnp from mnp using equation (8.2). 000 should be zero, so set 000 = 0.

e. Calculate the result vertex values of jkl by performing three sets of N^2 of size N Discrete

Fourier Series reversing step c.

You may wish to write 1-D and 2-Discrete Fourier Poisson solvers to get the details correct first.

Partial credit will be given for working 1-D and 2-D codes. Use the given 3-D charge distribution to

define the 2 D and 1-D charge distributions.

As the charge distribution is purely real and the resultant electric potential field is purely real, the

half of the calculation using imaginary values is wasted. It is possible to calculate the Discrete

Fourier Transform of Real Data efficiently when you are doing multiple transforms of the same size

by putting one set of real data into the real part of a size N complex vector and putting another set

of real data into the imaginary part. The Discrete Fourier Transform y of a purely real Data vector x

has the property: yN-k = yk* , the reversed second half of the transform array is the complex

conjugate of the first half. For a purely imaginary Data vector x, the y vector has the property:

yN-k = -yk* .

How long do steps c, d & e take for N = 16, 32, 48, 64, 80…. etc. ?

For what values of N does available computer memory become a problem?

In real modelling of atoms, molecules and semiconductors by this method , the distribution of

charges is averaged onto the nearest jkl grid vertex, the resultant electric potential at each grid

point jkl is calculated by the above algorithm and then each charge is moved by Maxwell’s laws in

the resultant Electric field. This process is repeated several thousand times for each model!

The results for the potential of our dipole brick may converge much faster if we set the faces to +

100, the edges to + 50 and the corners to + 25 (do you have any ideas why!). Test this for the 3-D

case and for any 2-D and 1-D cases if you have coded those up.

Plot j(N/2)(N/2), (N/2)k(N/2) & (N/2)(N/2)l for 0 < j,k,l < N on the same scale for your three largest

values of N. How long does this method of solving Poisson’s Equation take for N = (3,4,5) 24+k

0 < k < (as large as possible).? How would you expect the time taken to increase as N increases

(to leading order as a power of N!). For what values of N does available computer memory become

a problem?

Exercise #3 1 11 March 2013

Exercise #3 Part – II due 22:00pm Friday 3 May

(worth 15 % for M3SC and 40% for M 4/5 SC)

Background and discussion

The discrete Real Fourier series represents the periodic real data xj at N discrete points as the sum

of N cosine and sine functions in the form :

(II-1)

This would be a N/2-size discrete Complex Fourier Series if yk = ak + I bk . 0 < k < N/2

Although there are purely real algorithms to calculate xj given ak and bk and vice versa, most

commonly these problems are mapped into a procedure utilizing the Discrete Complex Fourier

Series and Transform.

One of these mappings take two real periodic data sets x1j and x2j 0 < j < N and packs them

into a periodic complex data set zj = x1j + x2j · I of size N and performs the size N Discrete

Complex Fourier Series on this data set: t = CN z . The resultant complex t vector can be

unpacked into the two real size N/2 Fourier Series coefficients by exploiting the fact that for a real

data series xj , the resultant complex yk have the property: yk = y*N-k . [II-2]

For a purely imaginary set of data, the resultant complex yk have the property: yk = - y*N-k .

Thus the yk for each real data set can be recovered for 0 < k < N/2 . In detail:

y 1k = ( tk + t*N-k )/N and y 2k = -I · ( tk - t*N-k )/N 0 < k < N/2.

To recover x1j and x2j 0 < j < N , given y1k and y2k 0 < k < N/2, extend y1k and y2k from

0 < k < N/2 to 0 < k < N using relation [II-2] and then create the complex vector

wk = y1k + y2k · I for 0 < k < N and perform the Discrete Fourier Transform CN-1 w .

The real and imaginary parts of the resultant complex array are the two desired real periodic data.

8. Write the C function to calculate the Discrete Fourier Series of two real data sets, each of

size N and find the two complex vectors z1 and z2 of size N/2 of their Fourier coefficients.

This function should have the prototype:

int TwoReal2Complex(double *R1, double *R2,

complex double *Z1,complex double *Z2,

complex double *w, complex double *Wp, int N)

where: *R1 is a pointer to the type double array to hold R1 from R1[0] to R1[N-1];

*R2 is a pointer to the type double array holding R2 from R2[0] to R2[N-1];

*Z1 is a pointer to the type complex double array holding Z1 from Z1[0] to Z1[N/2];

*Z2 is a pointer to the type complex double array holding Z2 from Z2[0] to Z2[N/2];

*w is a pointer to the type double array w[0] to w[your choice!] that can be used as

vector temporary storage for intermediate complex works required by your

FastCS(… function.

*Wp is a pointer to the vector or powers of WN created in the section 1 function above

N is the size of the desired transform.)

M 3/4/5 SC – SCIENTIFIC COMPUTING -- Due: 3 May 2013

Exercise #3 2 11 March 2013

9. Write the C function to calculate the Discrete Fourier Transform of two real periodic data sets,

each of size N given the two complex vectors z1 and z2 of size N/2 of their Fourier

coefficients This function should have the prototype:

int Complex2TwoReal(complex double *Z1,complex double *Z2,

double *R1, double *R2, complex double *w,

complex double *Wp, int N)

where

*Z1 is a pointer to the type complex double array holding Z1 from Z1[0] to Z1[N/2];

*Z2 is a pointer to the type complex double array holding Z2 from Z2[0] to Z2[N/2];

*R1 is a pointer to the type double array to hold R1 from R1[0] to R1[N-1];

*R2 is a pointer to the type double array holding R2 from R2[0] to R2[N-1];

*w is a pointer to the type double array w[0] to w[your choice!] that can be used as

vector temporary storage for intermediate complex works required by your

FastCS(… function.

*Wp is a pointer to the vector or powers of WN created in the section 1 function above

N is the size of the desired transform.)

10. Poisson’s Equation: 2 = - describes many fields in physics: electric field potential just

one example. A distribution of charge generates an electric field in which the charges

move in obeyance to Gauss’s laws (actually, Maxwell’s Laws) .

In Cartesian (x,y,z) co-ordinates this equation becomes

2

2

2

2

2

2

x

y

z

. (10.1)

If we write the vector r = x i + y j + z k and the vector K = kx i + ky j + kz k where i, j and k

are the usual unit vectors we can solve equation (8) by the separation of variables technique.

Let us assume that (x,y,z) = ( r) = exp( I K · r ) ( K) for some K vector. And assume that

(x,y,z) = ( r) = exp( I K · r ) ( K) . Substituting these expressions into

equation (8) leads to the relation

- (kx · kx + ky · ky + kz · kz) ( K) = - ( K) (10.2)

for every K vector.

Assume we want to solve equation (8.1) in a triply-periodic unit cube. This restricts kx, ky and kz

to take the discrete values: 2 π m, 2 π n and 2 π p for the integers m, n & p.

(We omit using o as an integer label as it may become confused with zero!)

Let (x,y,z) and (x,y,z) be defined at the vertices of a periodic cubical grid of size 1 with an

internal spacing of size = (1/N). So our functions and are defined at the points where

(x,y,z) take the values ( l/N, m/N, n/N), We can label each vertex value lmn and lmn.

This allows us to write : as a sum of periodic modes in each space dimension:

(x,y,z) = ( j, k , l ) =

1

, ,

0

N

mnp

m n p

exp(? ? ? ? ?

?

) exp(? ? ? ? ?

?

) exp(? ? ? ? ?

?

)

and (10.3)

(x,y,z) = ( j , k , l ) =

1

, ,

0

N

mnp

m n p

exp(? ? ? ? ?

?

) exp(? ? ? ? ?

?

) exp(? ? ? ? ?

?

)

for 0 < j, k, l < N .We assume a given periodic charge is assigned a value jkl at each vertex.

The problem is to solve Poisson’s equation (8.1) to find jkl at each (x,y,z) grid point. Both and

are assumed to have a zero mean value in our unit cube and to be real everywhere.

Find the electric potential ( l/N, m/N, n/N) for the charge distribution:

= -200 3/8 < x < 1/2 , 3/8 < y < 5/8, 7/16 < z < 9/16,

M 3/4/5 SC – SCIENTIFIC COMPUTING -- Due: 3 May 2013

Exercise #3 3 11 March 2013

= 0 x=0

= 200 1/2 < x < 5/8 , 3/8 < y < 5/8, 7/16 < z < 9/16,

You may assume that =0 everywhere else in the cube.

Find the solutions jkl for N = 16, 32, 48, 64, 80, 96,128, 160, 192, 256, 320 … until the problem

is too big for the PC.

Suggestions:

a. Allocate type double arrays of size N3 to hold the input, the output and the work arrays.

{I would recommend allocating these three arrays as single dimension arrays of size N3 .

Use the mapping p = j+ Nk + N2l to find the location of vertex j,k,l }

b. Set the jkl array to zero. For j from 3N/8 to N/2-1, for k from 3N/8 to 5N/8 and for

l from 7N/16 to 9N/16 set jkl = -200. Repeat for j from N/2+1 to 5N/8 with jkl = 200.

c. Find the Discrete Fourier transform of this distribution jkl by:

i. For k from 0 to N-1, For l from 0 to N-1, do the size N Discrete Fourier Transform

starting at 0kl with an initial skip of 1.

ii. For j from 0 to N-1, For l from 0 to N-1, do the size N Discrete Fourier Transform

starting at j0l with an initial skip of N.

iii. For j from 0 to N, For k from 0 to N, do the size N Discrete Fourier Transform

starting at jk0 with an initial skip of N2.

d. Calculate mnp from mnp using equation (8.2). 000 should be zero, so set 000 = 0.

e. Calculate the result vertex values of jkl by performing three sets of N^2 of size N Discrete

Fourier Series reversing step c.

You may wish to write 1-D and 2-Discrete Fourier Poisson solvers to get the details correct first.

Partial credit will be given for working 1-D and 2-D codes. Use the given 3-D charge distribution to

define the 2 D and 1-D charge distributions.

As the charge distribution is purely real and the resultant electric potential field is purely real, the

half of the calculation using imaginary values is wasted. It is possible to calculate the Discrete

Fourier Transform of Real Data efficiently when you are doing multiple transforms of the same size

by putting one set of real data into the real part of a size N complex vector and putting another set

of real data into the imaginary part. The Discrete Fourier Transform y of a purely real Data vector x

has the property: yN-k = yk* , the reversed second half of the transform array is the complex

conjugate of the first half. For a purely imaginary Data vector x, the y vector has the property:

yN-k = -yk* .

How long do steps c, d & e take for N = 16, 32, 48, 64, 80…. etc. ?

For what values of N does available computer memory become a problem?

In real modelling of atoms, molecules and semiconductors by this method , the distribution of

charges is averaged onto the nearest jkl grid vertex, the resultant electric potential at each grid

point jkl is calculated by the above algorithm and then each charge is moved by Maxwell’s laws in

the resultant Electric field. This process is repeated several thousand times for each model!

The results for the potential of our dipole brick may converge much faster if we set the faces to +

100, the edges to + 50 and the corners to + 25 (do you have any ideas why!). Test this for the 3-D

case and for any 2-D and 1-D cases if you have coded those up.

Plot j(N/2)(N/2), (N/2)k(N/2) & (N/2)(N/2)l for 0 < j,k,l < N on the same scale for your three largest

values of N. How long does this method of solving Poisson’s Equation take for N = (3,4,5) 24+k

0 < k < (as large as possible).? How would you expect the time taken to increase as N increases

(to leading order as a power of N!). For what values of N does available computer memory become

a problem?

下一篇：代写law paper