# Homework 4, due 23.59, 3/10-2016¶

## Curvilinear coordinates¶

The files relevant to this homework reside in `codes/fortran/HWK4`

.

In this homework you will compute derivatives and integrals on a logically rectangular domain (i.e. it has four sides and no holes) \(\Omega\). The strategy is to use a smooth mapping \((x,y)=(x(r,s),y(r,s))\) from the reference element \(\Omega_R = \{(r,s) \in [-1,1]^2\) } to \(\Omega\) to “transfer” the computations to the reference element.

For example if we want to approximate \(\frac{\partial u(x,y)}{\partial x}\) we use the chain rule to find

and similarly for

We can then introduce a Cartesian grid on the reference element:

```
hr = 2.d0/dble(nr)
hs = 2.d0/dble(ns)
do i = 0,nr
r(i) = -1.d0 + dble(i)*hr
end do
do i = 0,ns
s(i) = -1.d0 + dble(i)*hs
end do
```

to find \(u_r\) and \(u_s\) by standard finite difference formulas (see the subroutine `differentiate`

.) To find the metric \(r_x, r_y, s_x, s_y\) we can first compute \(x_r, x_s, y_r, y_s\) and then use the above formulas with \(u = x\) and \(u = y\). This yields

It is also possible to compute integrals on the reference element. For example, we have that

where \(J(r,s) = x_r y_s - x_s y_r\) is the surface element. The second integrand can be approximated on the reference element, for example with the trapezoidal rule.

Assignments:

Make and run the program

`homework4.f90`

and use the Matlab script (or some other plotting tool) to display the grid.Use the above formula to find (compute numerically) the metric \(r_x, r_y, s_x, s_y\).

Check that the results from 1) are correct by changing the mapping in

`xycoord.f90`

to map the reference element into some geometrical shape for which you know the area (for example a sector of an annulus) and compute it on the reference element using the trapezoidal rule.Compute \(u_x\) and \(u_y\) for some different functions \(u(x,y)\) and some different mappings. Write a subroutine that approximates the error

\begin{equation*} e_2(h_r,h_s) = \left(\int_{\Omega} \left(u_x(x,y)+u_y(x,y) - \left[(u_{\rm exact})_x + (u_{\rm exact})_y \right] \right)^2 dxdy \right)^{1/2}, \end{equation*}and plot the error as a function of the grid spacing for the different functions and the different mappings.

Use the chain rule to find an expression for \(\Delta u = u_{xx}+u_{yy}\). Discretize it and repeat the experiments above (optional).

Write up your findings neatly as a report and check it in to your repository.

## Sample results¶

Below we give some sample results, note that the errors are messured in the max-norm so the results you produce are not going to be identical. Also note that the errors are plotted as a function of an effective gridsize \(h_{\rm eff} = \sqrt{h_r h_s \max J}\) Here we use 3 combinations of grids and functions:

The three combinations are

```
! Combination 1
x_coord = r+0.1d0*s
y_coord = s
u = sin(xc)*cos(yc)
! Combination 2
x_coord = (2.d0+r+0.1d0*sin(5.d0*pi*s))*cos(0.5d0*pi*s)
y_coord = (2.d0+r+0.1d0*sin(5.d0*pi*s))*sin(0.5d0*pi*s)
u = exp(xc+yc)
! Combination 3
x_coord = r
y_coord = s + s*r**2
u = xc**2+yc**2
```

And the grids and results are: