<a href="https://colab.research.google.com/github/Jun-629/20MA573/blob/master/src/Hw9.ipynb" target="_parent"></a>

Consider 2-d PDE

\begin{align}
\frac{1}{2}\Delta v(x) - v(x) + x_1^2 + x_2^2 - x_1 - x_2 - \frac{3}{2} = 0, x \in O = (0,1)^2 \tag{1} \\
\end{align}

with its boundary data
$$v(x) = (x_1 - \frac{1}{2})^2 + (x_2 - \frac{1}{2})^2, x \notin O.$$
- Show that exact solution is
$$v(x) = (x_1 - \frac{1}{2})^2 + (x_2 - \frac{1}{2})^2.$$

__Pf:__

\begin{aligned}
L.H.S &= \frac{1}{2} \sum_{i=1}^2 \partial_{ii} v(x) - v(x) + x_1^2 + x_2^2 - x_1 - x_2 - \frac{3}{2} \\
&= \frac{1}{2} \cdot (2+2) - (x_1^2 - x_1 + \frac{1}{4}) - (x_2^2 - x_2 + \frac{1}{4}) + x_1^2 + x_2^2 - x_1 - x_2 - \frac{3}{2} \\
&= 1 - x_1^2 - x_2^2 + x_1 + x_2 - \frac{1}{2} + x_1^2 + x_2^2 - x_1 - x_2 - \frac{3}{2} \\
&= R.H.S \\
\end{aligned}

Also, $v(x)$ meets the boundary data, i.e., for $x \notin O$.

Therefore, the exact solution is 
$$v(x) = (x_1 - \frac{1}{2})^2 + (x_2 - \frac{1}{2})^2.$$

- Identify $\gamma, \ell^h, p^h$ in the CFD solution given by

\begin{align}
v(x) = \gamma \{ \ell^h(x) + \sum_{i=1}^d p^h(x+he_i|x)v(x+he_i) + p^h(x-he_i|x)v(x-he_i) \}. \tag{2} \\
\end{align}

__Soln:__

For CFD Solution, we already have that:

\begin{cases}
\partial_i v(x) = \frac{v(x+he_i) - v(x-he_i)}{2h} \\
\partial_{ii} v(x) = \frac{v(x+he_i) - 2v(x) + v(x-he_i)}{h^2} \\
\end{cases}

Then, denote $v = v(x)$, $v_i^+ = v(x+he_i)$ and $v_i^- = v(x-he_i)$. Take these equations into the equation $(1)$, we will have that:
$$\frac{1}{2} \sum_{i=1}^2 \frac{v_i^+ - 2v + v_i^-}{h^2} - v + x_1^2 + x_2^2 - x_1 - x_2 - \frac{3}{2} = 0 $$

\begin{aligned}
(\frac{2}{h^2} + 1)v &= \sum_{i=1}^2 \frac{1}{2h^2} v_i^+ + \sum_{i=1}^2 \frac{1}{2h^2} v_i^- + x_1^2 + x_2^2 - x_1 - x_2 - \frac{3}{2} \\
v(x) = v &= \frac{2}{2+h^2} \cdot \{ \frac{h^2}{2} (x_1^2 + x_2^2 - x_1 - x_2 - \frac{3}{2}) + \sum_{i=1}^2 \frac{1}{4} v_i^+ + \sum_{i=1}^2 \frac{1}{4} v_i^- \} \\
\end{aligned}

Thus, to compare the above equation to the equation $(2)$, we can deduce that

\begin{cases}
\gamma = \frac{2}{2+h^2}, \\
\ell^h(x) = \frac{h^2}{2}(x_1^2 + x_2^2 - x_1 - x_2 - \frac{3}{2}), \\
d = 2, \\
p^h(x \pm he_i|x) = \frac{1}{4}. \\
\end{cases}




- For $\{ h = 2^{-i}, i = 2,3,4,5 \}$, compute CFD solution and find maxnorm of error.


```
import numpy as np

def gamma(N):
  h = 1/N
  gamma = 2/(2 + h**2)
  return gamma

def l(N):
  h = 1/N
  l_h = np.zeros([N+1, N+1])
  for i in range(N+1):
    for j in range(N+1):
      if i*j == 0 or i == N or j == N:
        l_h[i,j] = h**2 / 2 * ((i/N - 0.5)**2 + (j/N - 0.5)**2)
      else:
        l_h[i,j] = h**2 / 2 * (((i+1)/N)**2 + ((j+1)/N)**2 - (i+1)/N - (j+1)/N - 3/2)
  return l_h

## define the iteration function F

def F(N, u):
  h = 1/N
  l_h = l(N)
  F = np.zeros([N+1, N+1])
  for i in range(0, N+1):
    for j in range(0, N+1):
      if i*j == 0 or i == N or j == N:
        F[i,j] = u[i,j]
      else:
        F[i,j] = gamma(N) * l_h[i,j] + gamma(N) * (u[i+1,j] + u[i,j+1] + u[i,j-1] + u[i-1,j]) / 4
  return F

# define the function of initial value of the PDE

def initial_value(N):
  Ini = np.zeros([N+1, N+1])
  for i in range(0, N+1):
    for j in range(0, N+1):
      if i*j == 0 or i == N or j == N:
        Ini[i,j] = ((i/N - 0.5)**2 + (j/N - 0.5)**2)
      else:
        Ini[i,j] = 0
  return Ini

def value_iteration(hat_eps, N):            # N: max iteration, hat_eps: tolerance
  v = initial_value(N)
  flag = 1
  n = 0                                     # the number of iteration

  eps = 1
  while flag:
    n += 1
    u = v
    v = F(N, u)
    eps = np.max(np.abs(u - v))
    if eps < hat_eps:
      flag = 0
  return [eps, n, v]

# Get the exact solution of the PDE

def exact_value(N):
  v_exact_value = np.zeros([N+1, N+1])
  for i in range(0, N+1):
    for j in range(0, N+1):
      v_exact_value[i, j] = (i/N - 0.5)**2 + (j/N - 0.5)**2
  return v_exact_value

v_exact_value = exact_value(4)

print("The exact value of the PDE is:")
print(v_exact_value)

# Find maxnorm of error

for i in range(2,6):
  error = np.max(np.abs(value_iteration(0.0001, 2**i)[2] - exact_value(2**i)))
  print("When h = ", 2**(-i), ", the maxnorn error is:")
  print(error)

```

    The exact value of the PDE is:
    [[0.5    0.3125 0.25   0.3125 0.5   ]
     [0.3125 0.125  0.0625 0.125  0.3125]
     [0.25   0.0625 0.     0.0625 0.25  ]
     [0.3125 0.125  0.0625 0.125  0.3125]
     [0.5    0.3125 0.25   0.3125 0.5   ]]
    When h =  0.25 , the maxnorn error is:
    0.02036086793376196
    When h =  0.125 , the maxnorn error is:
    0.007767404850966719
    When h =  0.0625 , the maxnorn error is:
    0.0048814207852426045
    When h =  0.03125 , the maxnorn error is:
    0.01847625984958414


c4d62ebf7a54aff3d700f3f3af590b19a94899da
ipynb
Jupyter Notebook
src/Hw9.ipynb
Jun-629/20MA573
addad663d2dede0422ae690e49b230815aea4c70
MIT
src/Hw9.ipynb
Jun-629/20MA573
addad663d2dede0422ae690e49b230815aea4c70
MIT
src/Hw9.ipynb
Jun-629/20MA573
addad663d2dede0422ae690e49b230815aea4c70
MIT
2020-02-05T21:42:08.000Z
2020-02-05T21:42:08.000Z
Qwen/Qwen-72B
1. YES
2. YES