这是一篇关于程序项目的数学计算代写，以下是作业具体要求：
Programming project in TMA4220 – part 1
The programming project will be split into two parts. The first part is an introduction into the finite element method and is designed to build up a solid code base for part 2, where you will actually solve larger, real life problems. For the second part you will have to choose one of several tasks to implement.
1 Gaussian quadrature
At the heart of every finite element code, lies the evaluation of an integral. This integral may be of varying complexity depending on the problem at hand, and many of these integrals does not even have a known analytical solution. Some integrals are possible to solve analytically, but of such computational complexity that it is impractical to do so. As such, one often refers to numerical integration schemes to do the core integration. One popular integration scheme is the Gaussian quadrature.
In one dimension the gauss quadrature takes the form
Z 1−1
g(z)dz ≈XNq
q=1
ρqg(zq),
where Nq is the number of integration points, zq are the Gaussian quadrature points and ρq are the associated Gaussian weights.
This extends to higher dimensions by
ZΩˆg(z)dz ≈XNq q=1ρqg(zq),and specifying the vector quadrature points zq as well as integrating over a suitable reference domain Ωˆ (i.e. squares or triangles in 2D, tetrahedra or cubes in 3D).
a) 1D quadrature
Write a Python function quadrature1D(a, b, Nq, g) that returns a value I where the variables are defined as:
I ∈ R value of the integral
a ∈ R integration start
b ∈ R integration end
Nq ∈ [1, 4] number of integration points
g : R → R function pointer∗
Verify that the function evaluates correctly by comparing with the analytical solution of the integral
Z 2
1exdx
Nq zq ρq
1-point-rule 0 2
2-point-rule −
pp
1/3 1
1/3 1
−p
3/5 5/9
3-point-rule 0 8/9
p
3/5 5/9
−q3+2√
6/57
18−√3036
4-point-rule −q3−2√6/5
718+√30q363−2√6/5
718+√30q363+2√6/5
718−√3036
Table 1: 1D gauss quadrature
b) 2D quadrature
Using all numerical quadratures, it is important to first map the function to the reference domain. In one dimension, this is the interval ζ ∈ [−1, 1]. In higher dimensions, we often map to barycentric coordinates (or area coordinates as they are known in 2D). The gauss points are then given as triplets in this coordinate system. The area coordinates are defined by
ζ1 =A1
Aζ2 =A2
Aζ3 =A3A
where A1, A2 and A3 are the area of the triangles depicted in figure 1 and A is the total area of the triangle. Note that these do not form a linear independent basis as ζ1 + ζ2 + ζ3 =1.
Nq (ζ1, ζ2, ζ3) ρ
1-point rule (1/3, 1/3, 1/3) 1
(1/2, 1/2, 0) 1/3
3-point rule (1/2, 0, 1/2) 1/3
(0, 1/2, 1/2) 1/3
(1/3, 1/3, 1/3) -9/16
4-point rule (3/5, 1/5, 1/5) 25/48
(1/5, 3/5, 1/5) 25/48
(1/5, 1/5, 3/5) 25/48
Table 2: 2D gauss quadrature
Write a Python function quadrature2D(p1, p2, p3, Nq, g) that returns a value I where
the variables are defined as:
Figure 1: Barycentric coordinates in two dimensions
I ∈ R value of the integral p1 ∈ R2
first corner point of the triangle p2 ∈ R2
second corner point of the triangle p3 ∈ R2
third corner point of the triangle Nq ∈ {1, 3, 4} number of integration points
g : R2 → R function pointer∗
Hint: An easy way of mapping barycentric coordinates ζ to physical coordinates x is by x =
ζ1p1 + ζ2p2 + ζ3p3, where pi, i = 1, 2, 3 is the corner points of the triangle.
Verify that the function evaluates correctly by comparing with the analytical solution of the integralZ ZΩlog(x + y) dx dy,where Ω is the triangle defined by the corner points (1, 0), (3, 1) and (3, 2).
(*) A function pointer in Python is a variable which represents a function instead of the usual
numerical values. In its simplest form it is declared as f = lambda x: x**2 + 1 which would cause the variable f to contain a pointer to the function f(x) = x 2 + 1. The function can then be evaluated by y = f(4) which should yield the result y = 17. A function may take in several arguments, i.e. f(x, y) = x2 + y 2 may be declared as f = lambda x, y: x**2 + y**2 again the evaluation of the function is straightforward y = f(2, 2)
Provided that the actual function body is capable of vector or matrix operations, then the input
arguments may be of vector or matrix form. The syntax remains unchanged by this.
2 Poisson in 2 dimensions
We are going to solve the two-dimensional Poisson problem, given by
∇2u(x, y) = −f(x, y) (1)
u(x, y)|r=1 = 0,
with f given in polar coordinates as
f(r, θ) = −8π cos
2πr2+ 16π2r2sin2πr2
on the domain Ω given by the unit disk, i.e. Ω =
(x, y) : x2 + y2 ≤ 1.
a) Analytical solution
Verify that the following expression is in fact a solution to the problem (3)
u(x, y) = sin2π(x2 + y2). (2)
b) Weak formulation
Show that the problem can be rewritten as
a(u, v) = l(v), ∀v ∈ X.
with the bilinear functional a and the linear functional l given by
a(u, v) = Z ZΩ
∇u · ∇v dx dy,
l(v) = Z ZΩ
fv dx dy.
What is the definition of the space X?
c) Galerkin projection
Instead of searching for a solution u in the entire space X we are going to be looking for a
solution in a much smaller space Xh ⊂ X. Let Ω be discretized into M triangles such that our
computational domain is the union of all of these Ω = ∪M
k=1Kk. Each triangle Kk is then defined
by its three corner nodes xi
. For each of these nodes there corresponds one basis function. The
space Xh is then defined by
Xh = {v ∈ X : v|Kk ∈ P1(Kk), 1 ≤ k ≤ M}
for which the basis functions {ϕi}
n
i=1 satisfy
Xh = span{ϕi}
n
i=1 ϕj (xi) =
and δij is the Kronecker delta. By searching for a solution uh ∈ Xh, it is then possible to write
this as a weighted sum of the basis functions, i.e. uh =
Pn
i=1 u
i
hϕi(x, y).
Show that the problem ”Find uh ∈ Xh such that a(uh, v) = l(v) ∀v ∈ Xh” is equivalent to the
following problem
Find u such that
Au = f (3)
with
A = [Aij ] = [a(ϕi
, ϕj )]
u = [u
i
h
]
f = [fi
] = [l(ϕi)].
d) Implementation
We are now going to actually solve the system (3). First we are going to take a look at the triangulation {Kk}. From the webpage https://wiki.math.ntnu.no/tma4220/2019h/
project you may download the mesh generators.
The function getDisk is generating the unit disk Ω. Plot at least three meshes of different sizes
using the mesh generated from this function.
e) Stiffness matrix
Build the stiffness matrix A. Use the Gaussian quadrature from exercise 1 to do this.
The matrix A should now be singular. Verify this in your code and explain why this is the case.
f) Right hand side
Build the right hand side vector f in the same manner as A.
g) Boundary conditions
Implement the homogeneous dirichlet boundary conditions. Describe what method you used for
this and how you did it.
h) Verification
Solve the system (3) and verify that you are getting (approximately) the same result as the analytical solution (2).
3 Neumann boundary conditions
We are going to change to boundary conditions of our problem to
∇2u(x, y) = −f(x, y) (4)
u(x, y)|∂ΩD
= 0,
∂u(x, y)
∂n
∂ΩN
= g(x, y),
with the source term f and exact solution u given as above, and g as
g(r, θ) = 4πr cos
2πr2
. (5)
The dirichlet boundary condition is defined on ∂ΩD =
x
2 + y
2 = 1, y < 0
, and the neumann
boundary condition as ∂ΩN =
x
2 + y
2 = 1, y > 0
shown in figure 2.
Figure 2: Dirichlet and Neumann boundary conditions
a) Boundary condition
Verify that (2) satisfies (5) at the boundary.
b) Variational formulation
How does a(·, ·) and l(·) change with the introduction of Neumann boundary conditions
c) Gauss quadrature
The neumann boundary condition is given as an integral and should be evaluated using Gaussian
quadrature. Modify your quadrature methods from task 1 to solve line integrals in two dimensions,
i.e. the method signature in quadrature1D(a, b, Nq, g) should change to a ∈ R
2
and b
∈ R
2
.
d) Implementation
Change your code from task 2 to solve this new boundary value problem. How does your solution
in the interior compare to the one you got in task 2? How does your solution at the boundary
compare?