# UC Berkeley Math 228B, Spring 2024, Problem Set 4

Due March 15

This content is protected and may not be shared, uploaded, or distributed.

1. Consider the boundary value problem \begin{aligned} u''''(x) &= f(x) \equiv 480x - 120, \qquad \text{for } x \in (0,1) &&(1)\\ u(0) &= u'(0) = u(1) = u'(1) = 0 &&(2) \end{aligned}

1. Derive the following Galerkin formulation for the problem (1)-(2) on some appropriate function space $$V_h$$: Find $$u_h \in V_h$$ such that \begin{aligned} \int_0^1 u_h''(x) v''(x)\, dx = \int_0^1 f(x) v(x)\, dx, \qquad \forall v \in V_h. &&(3) \end{aligned}

2. Define the triangulation $$T_h = \{K_1, K_2\}$$, where $$K_1 = [0,\frac12]$$ and $$K_2 = [\frac12,1]$$, and the function space \begin{aligned} V_h = \{ v\in C^1([0,1]) : v\big|_K \in \mathbb{P}_3(K)\ \forall K\in T_h,\ v(0)=v'(0)=v(1)=v'(1)=0 \}. \end{aligned} Find a basis $$\{\varphi_i\}$$ for $$V_h$$. Hint: Consider Hermite polynomials on each element.

3. Solve the Galerkin problem (3) using your basis functions. Plot the numerical solution $$u_h(x)$$ and the true solution $$u(x)$$.

2. Implement a Julia function with the syntax

   u = fempoi(p,t,e)

that solves Poissons’s equation $$-\nabla^2 u(x,y)=1$$ on the domain described by the unstructured triangular mesh p,t. The boundary conditions are homogeneous Neumann ($$n\cdot\nabla u=0$$) except for the nodes in the array e which are homogeneous Dirichlet ($$u=0$$).

Here are a few examples for testing the function:

# Square, Dirichlet left/bottom
pv = Float64[0 0; 1 0; 1 1; 0 1; 0 0]
p, t, e = pmesh(pv, 0.15, 0)
e = e[@. (p[e,1] < 1e-6) | (p[e,2] < 1e-6)]
u = fempoi(p, t, e)
tplot(p, t, u)

# Circle, all Dirichlet
n = 32; phi = 2pi*(0:n)/n
pv = [cos.(phi)  sin.(phi)]
p, t, e = pmesh(pv, 2pi/n, 0)
u = fempoi(p, t, e)
tplot(p, t, u)

# Generic polygon geometry, mixed Dirichlet/Neumann
x = 0:.1:1
y = 0.1*(-1).^(0:10)
pv = [x y; .5 .6; 0 .1]
p, t, e = pmesh(pv, 0.04, 0)
e = e[@. p[e,2] > (.6 - abs(p[e,1] - 0.5) - 1e-6)]
u = fempoi(p, t, e)
tplot(p, t, u)

3. Implement a Julia function with the syntax

   errors = poiconv(pv, hmax, nrefmax)

that solves the all-Dirichlet Poisson problem for the polygon pv, using the mesh parameters hmax and nref = 0,1,...,nrefmax. Consider the solution on the finest mesh the exact solution, and compute the max-norm of the errors at the nodes for all the other solutions (note that this is easy given how the meshes were refined – the common nodes appear first in each mesh). The output errors is a vector of length nrefmax containing all the errors.

Test the function using the commands below, which makes a convergence plot and estimates the rates:

hmax = 0.15
pv_square = Float64[0 0; 1 0; 1 1; 0 1; 0 0]
pv_polygon = Float64[0 0; 1 0; .5 .5; 1 1; 0 1; 0 0]

errors_square = poiconv(pv_square, hmax, 3)
errors_polygon = poiconv(pv_polygon, hmax, 3)
errors = [errors_square errors_polygon]

clf()
loglog(hmax ./ [1,2,4], errors)
rates = @. log2(errors[end-1,:]) - log2(errors[end,:])

Code Submission: Your Julia file needs to define the functions fempoi and poiconv, with exactly the requested names and input/output arguments, as well as any other supporting functions and variables that are required for your functions to run correctly.