[SOLVED] EC526 Assignment 3-Gaussian Elimination and Integration

Assignment 3: Gaussian Elimination and Integration

I Background

In the first part, you will write your own code to perform Gaussian elimination, extended to in- clude row pivoting to help with numerical stability. After a few warm up examples to test your code, you’ll use it to find the weights for arbitrary order Gaussian quadrature. This depends on the function int getLegendreZero(double* zero, double* a, int n) from the last problem set, if you couldn’t complete that problem please let us know! The zeros of the first 64 or so Legendre polynomials can be found on the web, so if need be you can proceed using those in the short run. https://pomax.github.io/bezierinfo/legendre-gauss.html

In the second part, you’ll put all of the pieces together and write your own numerical integration code. You’ll integrate a few functions we supply on the interval [−1, 1] comparing Gaussian integration with N terms in the sum to the trapezoid rule from Assignment 3 with h = 1/(2N).

In the third part you are going to push to multidimensional integrals and one 4d integral to run in parallel using OpenACC.

I.1 Warm up on Gaussian Elimination.

Gaussian elimination is a general algorithm to systematically solve systems of equations. In the lecture notes, we gave a blow by blow example for 3 unknowns:

−x + 2y − 5z = 17
2x + y + 3z = 0
4x − 3y + z = −10 (1)

This form lends itself nicely to multi-dimensional arrays. The left hand side can be stored as a 3×3 two-dimensional array, and the right hand side can be stored as a one-dimensional array of length 3. These are commonly referred to as A and b, respectively. These operations can all be done in place in A and b (with a temporary variable for interchanging rows when you pivot, and with the assumption that you’re not worried about saving the original matrix A). When the algorithm is done, b contains the solution!

Your task for this problem to solve a system of linear equations using Gaussian elimination. The routine needs three inputs:

1

• A 2D array A,
• A 1D array b,
• the dimension of the matrix and the array, N.

As noted above, your algorithm should be done in place, destroying the original contents of A and b. At the end of the algorithm, b will contain the solution of the linear system. A benefit of an in place algorithm is that you avoid allocating additional memory and carrying around additional references. While this isn’t an issue with problems this small, in Big Data applications spurious copies of data can be a killer. In this small data example you many want to save one copy of A and b outside the function to have the original data for future reference and testing. The function you need to program is:

where

• N is the dimension of the problem, that is, the number of equations and variables.
• A is a N by N 2D array which contains the coefficients of a system of equations.
• b is a 1D array of length N which contains the right hand side of a system of equations.

For example, for the system of equations in Eq. 1, one might write the code to form A and b, call the gaussian elimination function, and print the results as:

int gaussianElimination(double** A, double* b, int N);

which should print out

1.0 4.0 -2.0

2

Another fun example to also test doubles precision arithmetic. Mathematica can give the exact solution to

as rationals

x/16 + 3y + z = 2 x/5 + y + 2z = 1/4

5x − 3y + z = 13 (2) {x− > 1000/317, y− > 1011/1268, z− > −(747/1268)} (3)

Why is the solution 2 rational. By the way it is not too hard introduce a rational class (e.g. data type) on C++ do this. Fun but of course Mathematica is all set with this (see code at https://github.com/brower/EC526 2022/HW3code)

II Coding Exercise #1: Applied to Uniform Grid
Ok suppose we are not as smart as Gauss – who is!. (This is a warm up for the Gaussian integral

next.) Lets to the integral

x’s. are

or the equation to solve

􏰀N 􏰆1 (xj)i−1wj =

(xj)i = (1 − (−1)i)/i) (7)

j=1

N

−1

􏰆1 􏰀N f(x)dx ≃

wjf(xj) (4) on a uniform grid of N-points (or N − 1 intervals since I am assuming that end points x1 = −1 and

xN =1areincluded! ) ofwidthh=2/(N−1)
xj =−1,−1+h,−1+2h,···1−(N−2)h,1 (5)

That is the grid pints are
xj =−1+2(j−1)/(N−1) for j=1,2,3,···N (6)

SincethepositionsofthepointsarefixedwecanonlydemandexactanswersforNterms: 1,x,x2,···xN−1. JUST the same algebrix step of full Gaussian integration except you don’t have to precompute the

􏰀A[i][j]w[j] = b[i] j=1

where

i = 1,··· ,N (8)

−1

j=1

3

where have defined the matrix and vector
A[i][j] = (−1 + 2(j − 1)/(N − 1))i−1 and b[i] = (1 − (−1)i)/i) (9)

SEE mathematica notebook on GitHUB for the exact values of the weights!

It turns out the even N is easiery to do beause there is no point at x = 0 and Mathematica can’t figure out how to compute x0 = 1 for x = 0. The general N case neeeds to add this the first row b hand. Then n = 1 it is the Reimann itegral rule. For n = 2 the trapazoid rule for n = 3 simplson rule. But you can keep inventing higher order rules for any N. Pretty good even if you are not as smart as Gauss.

III Coding Exercise #2: Applied to Gaussian Quadrature

You should refer to our lecture notes on numerical integration for our base discussion of approximating integrals from −1 to 1 in the following form:

where have defined the matrix and vector A[i][j] = xi−1 and

b[i] =

􏰆1 −1

xi−1dx = (1 − (−1)i)/i (16)

for i,j = 1,2,···N.

􏰆1 􏰀N f(x)dx ≃

wif(xi) (10)

−1

i=1

We learned that we needed N points to integrate a polynomial of degree x2N−1 exactly. More im- portantly, we discussed how the points xi, i = 1,2,…,N exactly coincided with the zeroes of the Legendre polynomial PN (x), which we found in Problems Set #3.

Here we find the wights wi for each zero using Gaussian elimination to do Gaussian quadrature. (Yes! Everything seems to have Gauss’ name on it! He even invented the FFT.) The weights must obey the conditions,

􏰆1 −1

􏰆1 −1

x0 dx (11) x1 dx (12)

(13) xN−1dx (14)

−1 Sothattheweightsw1,w2,···wN arefoundbysolvingtheNlinearequations:

w1 +w2 +···+wN = x1w1 +x2w2 +···+xNwN =

. xN−1w +xN−1w +···+xN−1w =

􏰆1 1122 NN

N

􏰀A[i][j]w[j] = b[i] j=1

where i = 1,···N (15)

j

4

You will write a program gauss quad weight.c which, given a value N that you ask the user for, prints out the values of the zeroes xi and the weights wi. You should reuse the program getZeros.c from the previous assignment to find the zeroes, xi, of the Legendre polynomial PN (x). You can then set up the system of equations above and use the function gaussianElimination that you wrote in the previous part of the problem to compute the weights, wi.

As an example of how the code may work, let’s consider asking for the zeroes and weights for N = 3. User input is given on lines beginning with >.

You can assume we will never ask for N larger than 32. You should take it upon yourself to check your answers against Wikipedia: https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Legendre_quadrature.

III.1 Exercise #4. Comparing Uniform vs to Gausssian Grid Integrals

It’s now time to put the pieces together and integrate a few functions on the interval [−1, 1]. This depends heavily on the code you wrote in the previous problem. Modify the integration program to again evaluate

􏰆1
x8 dx =?

−1 􏰆1

cos(πx/2) dx =? −1

􏰆1 1 dx=? (17) −1 x2 + 1

5

but this time use Gaussian integration for N = 2, 4, 8. · · · , s32 and plot error in comparison with the trapezoidal rule. You will want to keep the Trapezoid rule method to be able to make comparisons between Gaussian integration for N point with Trapezoidal rule for 2N points.

Next extend this to multi-dimensional integrals (The next problem we will do some really hard multiple integral using GPUs and parallelization with OpenACC!)

As a quick refresher, a multiple dimensional integrals

􏰆􏰆 􏰆1􏰆1
g(x, y)dA = dy dx g(x, y) (18)

􏰈 −1−1
Applying Gaussian integration to each integral is straight forward,

􏰆1 􏰆1 􏰆1 􏰀N dy dx g(x, y) ≈ dy

i=1

−1−1 −1 NN

wig(xi, y) (19) ≈ 􏰀􏰀wiwjg(xi,yj) (20)

j=1 i=1

This is just a nested loop for a 1d integrals over the x-coordinate and the y-coordinated. Note that the weights wi and the value of xi and yj are exactly the same as in the one dimensional example. This approximation generalizes trivially and dimension on a “hyper-cubic” cell.

III.2 Integrating a few multiple double Integrals

This problem should seem familiar: we’re going to integrate a few problems on the interval [−1, 1], ex- cept in multiple dimensions. Write a main program test integrate 2d.c that performs the following

6

three integrals using Gaussian integration for N = 2 to 24 in the sum:

􏰆1 􏰆1
dy dx[x8 +y8 +(y−1)3(x−3)5]=? (21)

−1 −1

􏰆1 􏰆1

dy dx 􏰂x2 − y2 + 2 =? (22) −1

−1
􏰆1􏰆12y2 π

dy dxe−x−8 cos(πx)sin( y)=? (23)

8
You can test you code against Mathematica. (Don’t be shy about getting help with Mathematica!)

−1 −1

7