\documentclass[12pt]{article}
\usepackage{amsmath, amssymb, amsthm}
\usepackage{indentfirst} 
\usepackage{tabularx}
\usepackage{graphicx} 
\usepackage{float}
\usepackage{esint}
\usepackage[left=1in, right=1in]{geometry}
\usepackage[T1]{fontenc}
\usepackage{amsmath} % used for boldsymbol.
\geometry{}
\usepackage{filecontents}
\begin{filecontents*}{harmonicfunction.bib}
@book{strauss_partial_2007,
	title = {Partial {Differential} {Equations}: {An} {Introduction}},
	isbn = {978-0-470-05456-7},
	shorttitle = {Partial {Differential} {Equations}},
	abstract = {Partial Differential Equations presents a balanced and comprehensive introduction to the concepts and techniques required to solve problems containing unknown functions of multiple variables. While focusing on the three most classical partial differential equations (PDEs)—the wave, heat, and Laplace equations—this detailed text also presents a broad practical perspective that merges mathematical concepts with real-world application in diverse areas including molecular structure, photon and electron interactions, radiation of electromagnetic waves, vibrations of a solid, and many more.Rigorous pedagogical tools aid in student comprehension; advanced topics are introduced frequently, with minimal technical jargon, and a wealth of exercises reinforce vital skills and invite additional self-study. Topics are presented in a logical progression, with major concepts such as wave propagation, heat and diffusion, electrostatics, and quantum mechanics placed in contexts familiar to students of various fields in science and engineering. By understanding the properties and applications of PDEs, students will be equipped to better analyze and interpret central processes of the natural world.},
	language = {en},
	publisher = {John Wiley \& Sons},
	author = {Strauss, Walter A.},
	month = dec,
	year = {2007},
	note = {Google-Books-ID: m2hvDwAAQBAJ},
	keywords = {Mathematics / Differential Equations / Partial, Mathematics / Mathematical Analysis, thesis}
}

@book{evans_partial_2010,
	title = {Partial {Differential} {Equations}},
	isbn = {978-0-8218-4974-3},
	abstract = {This is the second edition of the now definitive text on partial differential equations (PDE). It offers a comprehensive survey of modern techniques in the theoretical study of PDE with particular emphasis on nonlinear equations. Its wide scope and clear exposition make it a great text for a graduate course in PDE. For this edition, the author has made numerous changes, including a new chapter on nonlinear wave equations, more than 80 new exercises, several new sections, a significantly expanded bibliography. About the First Edition: I have used this book for both regular PDE and topics courses. It has a wonderful combination of insight and technical detail. ... Evans' book is evidence of his mastering of the field and the clarity of presentation. --Luis Caffarelli, University of Texas It is fun to teach from Evans' book. It explains many of the essential ideas and techniques of partial differential equations ... Every graduate student in analysis should read it. --David Jerison, MIT I usePartial Differential Equationsto prepare my students for their Topic exam, which is a requirement before starting working on their dissertation. The book provides an excellent account of PDE's ... I am very happy with the preparation it provides my students. --Carlos Kenig, University of Chicago Evans' book has already attained the status of a classic. It is a clear choice for students just learning the subject, as well as for experts who wish to broaden their knowledge ... An outstanding reference for many aspects of the field. --Rafe Mazzeo, Stanford University},
	language = {en},
	publisher = {American Mathematical Soc.},
	author = {Evans, Lawrence C.},
	year = {2010},
	note = {Google-Books-ID: Xnu0o\_EJrCQC},
	keywords = {Mathematics / Differential Equations / General, thesis}
}

@book{rudin_principles_1976,
	title = {Principles of {Mathematical} {Analysis}},
	isbn = {978-0-07-085613-4},
	abstract = {The third edition of this well known text continues to provide a solid foundation in mathematical analysis for undergraduate and first-year graduate students. The text begins with a discussion of the real number system as a complete ordered field. (Dedekind's construction is now treated in an appendix to Chapter I.) The topological background needed for the development of convergence, continuity, differentiation and integration is provided in Chapter 2. There is a new section on the gamma function, and many new and interesting exercises are included. This text is part of the Walter Rudin Student Series in Advanced Mathematics.},
	language = {en},
	publisher = {McGraw-Hill},
	author = {Rudin, Walter},
	year = {1976},
	note = {Google-Books-ID: kwqzPAAACAAJ},
	keywords = {Mathematics / Calculus, thesis}
@misc{jiang_roubles_nodate,
	title = {Rouble's woes reap mixed harvest[1]- {Chinadaily}.com.cn},
	url = {http://www.chinadaily.com.cn/business/2015-04/13/content_20416109.htm},
	urldate = {2019-04-08},
	author = {Jiang, XueQing},
	file = {Rouble's woes reap mixed harvest[1]- Chinadaily.com.cn:C\:\\Users\\Yilin Jiang\\Zotero\\storage\\A6IELEJC\\content_20416109.html:text/html}
}
@misc{gungahlin_weather_centre_weather_nodate,
	title = {Weather station to be built on {Mt}. {Everest}!},
	url = {http://www.gungahlinweather.com/articles/everest.htm},
	urldate = {2019-04-24},
	author = {Gungahlin Weather Centre},
	file = {Weather station to be built on Mt. Everest!:C\:\\Users\\Yilin Jiang\\Zotero\\storage\\CGHT3Z2W\\everest.html:text/html}
}
@inproceedings{shepard_two-dimensional_1968,
	address = {Not Known},
	title = {A two-dimensional interpolation function for irregularly-spaced data},
	url = {http://portal.acm.org/citation.cfm?doid=800186.810616},
	doi = {10.1145/800186.810616},
	language = {en},
	urldate = {2019-04-24},
	booktitle = {Proceedings of the 1968 23rd {ACM} national conference on   -},
	publisher = {ACM Press},
	author = {Shepard, Donald},
	year = {1968},
	keywords = {thesis},
	pages = {517--524},
	file = {Shepard - 1968 - A two-dimensional interpolation function for irreg.pdf:C\:\\Users\\Yilin Jiang\\Zotero\\storage\\ENJ2D327\\Shepard - 1968 - A two-dimensional interpolation function for irreg.pdf:application/pdf}
}
@article{meijering_chronology_2002,
	title = {A chronology of interpolation: from ancient astronomy to modern signal and image processing},
	volume = {90},
	issn = {00189219},
	shorttitle = {A chronology of interpolation},
	url = {http://ieeexplore.ieee.org/document/993400/},
	doi = {10.1109/5.993400},
	abstract = {This paper presents a chronological overview of the developments in interpolation theory, from the earliest times to the present date. It brings out the connections between the results obtained in diﬀerent ages, thereby putting the techniques currently used in signal and image processing into historical perspective. A summary of the insights and recommendations that follow from relatively recent theoretical as well as experimental studies concludes the presentation.},
	language = {en},
	number = {3},
	urldate = {2019-04-08},
	journal = {Proceedings of the IEEE},
	author = {Meijering, E.},
	month = mar,
	year = {2002},
	pages = {319--342},
	file = {Meijering - 2002 - A chronology of interpolation from ancient astron.pdf:C\:\\Users\\Yilin Jiang\\Zotero\\storage\\39AEEX6W\\Meijering - 2002 - A chronology of interpolation from ancient astron.pdf:application/pdf}
}
@article{maude_interpolationmainly_1973,
	title = {Interpolation—mainly for graph plotters},
	volume = {16},
	issn = {0010-4620},
	url = {https://academic.oup.com/comjnl/article/16/1/64/435132},
	doi = {10.1093/comjnl/16.1.64},
	abstract = {Abstract.  A method of interpolation is given in which the function and its first derivative are continuous at all points, but which is simpler than spline inte},
	language = {en},
	number = {1},
	urldate = {2019-04-24},
	journal = {The Computer Journal},
	author = {Maude, A. D.},
	month = jan,
	year = {1973},
	keywords = {thesis},
	pages = {64--65},
	file = {Full Text PDF:C\:\\Users\\Yilin Jiang\\Zotero\\storage\\2HKAJU8V\\Maude - 1973 - Interpolation—mainly for graph plotters.pdf:application/pdf;Snapshot:C\:\\Users\\Yilin Jiang\\Zotero\\storage\\NM26FJDD\\435132.html:text/html}
}
@article{birkhoff_interpolation_1973,
	title = {Interpolation to boundary data in triangles},
	volume = {42},
	issn = {0022-247X},
	url = {http://www.sciencedirect.com/science/article/pii/0022247X73901546},
	doi = {10.1016/0022-247X(73)90154-6},
	number = {2},
	urldate = {2019-04-22},
	journal = {Journal of Mathematical Analysis and Applications},
	author = {Birkhoff, Garrett},
	month = may,
	year = {1973},
	pages = {474--484},
	file = {ScienceDirect Full Text PDF:C\:\\Users\\Yilin Jiang\\Zotero\\storage\\56Y4TXYL\\Birkhoff - 1973 - Interpolation to boundary data in triangles.pdf:application/pdf;ScienceDirect Snapshot:C\:\\Users\\Yilin Jiang\\Zotero\\storage\\2XIER6VY\\0022247X73901546.html:text/html}
}

@article{franke_scattered_1982,
	title = {Scattered data interpolation: tests of some methods},
	volume = {38},
	issn = {0025-5718},
	shorttitle = {Scattered data interpolation},
	url = {http://www.ams.org/jourcgi/jour-getitem?pii=S0025-5718-1982-0637296-4},
	doi = {10.1090/S0025-5718-1982-0637296-4},
	abstract = {This paper is concerned with the evaluation of methods for scattered data interpolation and some of the results of the tests when applied to a number of methods. The process involves evaluation of the methods in terms of timing, storage, accuracy, visual pleasantness of the surface, and ease of implementation. To indicate the flavor of the type of results obtained, we give a summary table and representative perspective plots of several surfaces.},
	language = {en},
	number = {157},
	urldate = {2019-04-23},
	journal = {Mathematics of Computation},
	author = {Franke, Richard},
	month = jan,
	year = {1982},
	keywords = {thesis},
	pages = {181--181},
	file = {Franke - 1982 - Scattered data interpolation tests of some method.pdf:C\:\\Users\\Yilin Jiang\\Zotero\\storage\\DWKRNUCF\\Franke - 1982 - Scattered data interpolation tests of some method.pdf:application/pdf}
}

@article{mclain_two_1976,
	title = {Two {Dimensional} {Interpolation} from {Random} {Data}},
	volume = {19},
	issn = {0010-4620, 1460-2067},
	url = {https://academic.oup.com/comjnl/article-lookup/doi/10.1093/comjnl/19.2.178},
	doi = {10.1093/comjnl/19.2.178},
	language = {en},
	number = {2},
	urldate = {2019-04-23},
	journal = {The Computer Journal},
	author = {McLain, D. H.},
	month = may,
	year = {1976},
	keywords = {thesis},
	pages = {178--181},
	file = {McLain - 1976 - Two Dimensional Interpolation from Random Data.pdf:C\:\\Users\\Yilin Jiang\\Zotero\\storage\\UU84RA8J\\McLain - 1976 - Two Dimensional Interpolation from Random Data.pdf:application/pdf}
}

@article{ohashi_astronomy_nodate,
	title = {{ASTRONOMY} {AND} {MATHEMATICS} {OF} {YIXING}},
	language = {en},
	author = {Ôhashi, Yukio},
	keywords = {thsis},
	pages = {8},
	file = {Ôhashi - ASTRONOMY AND MATHEMATICS OF YIXING.pdf:C\:\\Users\\Yilin Jiang\\Zotero\\storage\\WC8S8ZAQ\\Ôhashi - ASTRONOMY AND MATHEMATICS OF YIXING.pdf:application/pdf}
}

@article{wong_rainfall_2003,
	title = {Rainfall prediction model using soft computing technique},
	volume = {7},
	issn = {1432-7643, 1433-7479},
	url = {http://link.springer.com/10.1007/s00500-002-0232-4},
	doi = {10.1007/s00500-002-0232-4},
	language = {en},
	number = {6},
	urldate = {2019-04-23},
	journal = {Soft Computing - A Fusion of Foundations, Methodologies and Applications},
	author = {Wong, K. W. and Wong, P. M. and Gedeon, T. D. and Fung, C. C.},
	month = may,
	year = {2003},
	pages = {434--438},
	file = {Wong 等。 - 2003 - Rainfall prediction model using soft computing tec.pdf:C\:\\Users\\Yilin Jiang\\Zotero\\storage\\DDU9CFDG\\Wong 等。 - 2003 - Rainfall prediction model using soft computing tec.pdf:application/pdf}
	}
@book{sauer_numerical_2011,
	address = {USA},
	edition = {2nd},
	title = {Numerical {Analysis}},
	isbn = {978-0-321-78367-7},
	abstract = {Normal 0 false false false Numerical Analysis, Second Edition, is a modern and readable text. This book covers not only the standard topics but also some more advanced numerical methods being used by computational scientists and engineerstopics such as compression, forward and backward error analysis, and iterative methods of solving equationsall while maintaining a level of discussion appropriate for undergraduates. Each chapter contains a Reality Check, which is an extended exploration of relevant application areas that can launch individual or team projects. MATLAB is used throughout to demonstrate and implement numerical methods. The Second Edition features many noteworthy improvements based on feedback from users, such as new coverage of Cholesky factorization, GMRES methods, and nonlinear PDEs.},
	publisher = {Addison-Wesley Publishing Company},
	author = {Sauer, Timothy},
	year = {2011},
	keywords = {thesis}
}

@book{sauter_boundary_2011,
	address = {Berlin Heidelberg},
	series = {Springer {Series} in {Computational} {Mathematics}},
	title = {Boundary {Element} {Methods}},
	isbn = {978-3-540-68092-5},
	url = {https://www.springer.com/us/book/9783540680925},
	abstract = {This work presents a thorough treatment of boundary element methods (BEM) for solving strongly elliptic boundary integral equations obtained from boundary reduction of elliptic boundary value problems in IR3. The book is self-contained, the prerequisites on elliptic partial differential and integral equations being presented in Chapters 2 and 3. The main focus is on the development, analysis, and implementation of Galerkin boundary element methods, which is one of the most flexible and robust numerical discretization methods for integral equations. For the efficient realization of the Galerkin BEM, it is essential to replace time-consuming steps in the numerical solution process with fast algorithms. In Chapters 5-9 these methods are developed, analyzed, and formulated in an algorithmic way.},
	language = {en},
	urldate = {2019-05-09},
	publisher = {Springer-Verlag},
	author = {Sauter, Stefan and Schwab, Christoph},
	year = {2011},
	keywords = {thesis},
	file = {Snapshot:C\:\\Users\\Yilin Jiang\\Zotero\\storage\\WBNHNLHX\\9783540680925.html:text/html}
}

\end{filecontents*}
\usepackage[numbers]{natbib}
\usepackage{url}

\newtheorem{thm}{Theorem}[section]
\newtheorem{Def}{Definition}[section]
\newtheorem{lemma}{Lemma}
\renewcommand{\proofname}{\textbf{Proof}}
\renewcommand{\vec}[1]{\boldsymbol{#1}} % Uncomment for BOLD vectors.

\newcommand{\eat}[1]{}

\title{Partial Differential Equations in Data Analysis}
\author{Yilin Jiang}
\date{ }




\begin{document}

\begin{titlepage}
   \begin{center}
       \vspace*{1cm}
 
       UNIVERSITY OF OKLAHOMA
 
       \vspace{0.3cm}
       GRADUATE COLLEGE
 
       \vspace{2.5cm}
 
       Partial Differential Equations in Data Analysis
 
       \vfill
 
       A THESIS\\
       SUBMITTED TO THE GRADUATE FACULTY\\
       \vspace{0.3cm}
       in partial fulfilment of the requirements for the 
       \vspace{0.3cm}
       Degree of
 
       \vspace{0.3cm}
       MASTER OF SCIENCE \\
       \vspace{1.5cm}
       BY\\
       \vspace{0.3cm}
       Yilin Jiang\\
       Norman, Oklahoma\\
       2019
    \end{center}
\end{titlepage}
\begin{titlepage}
   \begin{center}
       \vspace*{1cm}
 
        Partial Differential Equations in Data Analysis
 
       \vspace{1cm}
       A THESIS APPROVED FOR THE \\
       DEPARTMENT OF MATHEMATICS 
 
       \vfill
 
       BY\\
       \end{center}
\begin{flushright}      
Dr. Meijun Zhu, Chair\\
\vspace{0.2cm}
Dr. John Albert\\
\vspace{0.2cm}
Dr. Pengfei Zhang
\end{flushright}
\end{titlepage}


\clearpage
\thispagestyle{empty}
\begin{center}
\vspace*{\fill}

\textcopyright\ Copyright by Yilin Jiang 2019\\
 All rights reserved.
\end{center}
\clearpage

\pagenumbering{roman}
\setcounter{page}{4}
\tableofcontents
\newpage
\thispagestyle{plain}
\begin{center}
    \Large
    \textbf{Partial Differential Equations in Data Analysis}
 
    
    \vspace{0.4cm}
    \textbf{Yilin Jiang}
 
 	\vspace{0.4cm}
     \begin{abstract}
In this thesis, we introduce new interpolation methods for  two dimensional data via constructing harmonic functions passing through the given data. Two ways to construct the harmonic function are introduced: (1) constructing the harmonic function via the heat equation,  and (2) constructing the harmonic function via the boundary element method 
\end{abstract}

\end{center}
\newpage
\eat{\maketitle}
\pagenumbering{arabic}
\setcounter{page}{1}

\section{Introduction}
Weather forecasting is not an easy job, and sometimes its accuracy can be affected by many factors. In middle of 2014, with the sharp fall of oil prices, the Russian rouble had declined drastically. However, at beginning of 2015, China’s National Weather Center said in its own microblog that China’s weather forecasts have become less accurate due to the Russian financial crisis \cite{jiang_roubles_nodate}.  Because of the funding cut, Russian’s weather center reduced its upper air data collecting through radiosonde observations from twice a day to once a day. Because of geographical reason, accurate upper air weather data from Russia are crucial to China’s weather forecast especially for predicting cold air from Siberia. Missing data (or in this case, simply less data) problems are not rare in real life data analysis. Interpolation methods are some of the most useful tools to deal with such a problem.The main idea of interpolation methods is to find an approximating function that passes through all the given data points. 
\begin{Def}
A function $P:\mathbb{R}^n \longrightarrow \mathbb{R} $ satisfying
\begin{align}
P(x_i)\, = \,y_i\quad  i = 1,2,3,...,n,
\end{align}
is called an interpolating function passing through distinct $(x_1,y_1),(x_2,y_2),...,(x_n,y_n)$, where $x_1,x_2,...x_n \in \mathbb{R}^2, y_1,y_2,...y_n \in \mathbb{R}$.  $(x_1,y_1),(x_2,y_2),...,(x_n,y_n)$ are called interpolating nodes.
\end{Def}
Using interpolation methods, one can fill in missing data on a given interval by calculating the values of the interpolating function. Hence, a more accurate weather forecast can be expected. Another scenario is that there are many tough places to forecast weather around the world. One may build a personal weather station at home easily, but building weather station in some places is a time and money consuming task. The first weather station on Mountain Everest was not built untill 2004, 50 years after the first climbers stood at the top of the mountain, in 1953. \cite{gungahlin_weather_centre_weather_nodate}. In these cases, data are only given on the boundary, and the interpolation methods are needed to estimate the data inside the given area. 

Before some well known interpolation methods like Lagrange method and Newton's method were invented during the Age of Scientific Revolution, many interpolation methods were used in ancient China. Second order interpolation was pioneered by the Chinese early astronomer Liu Zhuo who used this technique to calculate the relative distance between the sun and the moon \citep{meijering_chronology_2002}. However, the mathematician, astronomer, and monk Yixing created a better interpolation method based on Liu Zhuo's method, called Da Yan Calendar, using unequal interval second order interpolation \cite{ohashi_astronomy_nodate}. In fact, many mathematicians tried to find a reasonable way to find an approximation function interpolating the given randomly distributed data on a two dimensional plane. Subdividing the domain into small regions like triangles or rectangles  and connecting the given data points is a common interpolation method \citep{birkhoff_interpolation_1973} \citep{maude_interpolationmainly_1973} \citep{mclain_two_1976}. In \cite{maude_interpolationmainly_1973}, Maude introduced a way that divided the domain into several regions formed by intersections of circles. In \cite{mclain_two_1976}, however, the author revealed that Maude's result will be awkward when some special number of data points are given. In \cite{shepard_two-dimensional_1968}, Shepard came up with an algorithm called the Inverse Distance Weight interpolation. In his algorithm, the interpolation function had the form $f(x) = \Sigma_{i=1}^{n} w_i(x)f_i$, where $w_i(x)$ is the weight coefficient which depends on the reciprocal of distance between interpolate point $x$ and all other $n$ known points $x_i$. A group of scientists studied Shepard's method and improved the interpolation function through modifying the weight function $w_i$. However, in \citep{franke_scattered_1982}, Franke pointed out that Shepard's method is very dependent on weight function. Franke then listed several unacceptable weight functions and time consuming functions. Both the subdividing method and the Inverse Distance Weight method require much work to be done before interpolating. Our method focus on providing another way to interpolate the given data points. 

\section{Interpolation} 
There are many ways to find $P$. Common interpolation methods include Lagrange  Interpolation, Newton Interpolation, Hermite Interpolation, Linear Interpolation, and Cubic Spline Interpolation. Five interpolation methods are shown in Figure 1, where they are all used to find the approximation function of $f(x) = \frac{1}{x^2+1},x \in [-5,5]$. Notice that the node points in Figure 1 are equidistant nodes.
\subsection{One Dimensional Interpolation}
\begin{enumerate}
\item
Polynomial Interpolation\\
In 1885, a very famous theorem proved by Karl Weierstrass indicated the use of polynomial to approximate a continuous  function:
\begin{thm}\textup{(Weierstrass Approximation Theorem)}\\
Suppose $f$ is a continuous real-valued function on $ [a,b] $. For any given $\epsilon > 0$, there exists a polynomial $P$ on $ [a,b] $ such that
$|f(x) - P(x)| < \epsilon $ for all $x \in [a,b]$. 
\end{thm}
\setlength{\parindent}{1em} In words, any continuous function on an interval can be uniformly approximated by the polynomial to any degree of accuracy.\\
If $P(x)$ satisfies (1) and is a polynomial degree less than $n$, i.e. 
\[
P_n(x) = a_0+a_1x+a_2x^2+...+a_nx^n,
\]
where $a_i \in \mathbb{R}$,\,then $P(x)$ is the interpolating polynomial, the corresponding method is called the polynomial interpolation. The following theorem shows that such a polynomial exists and it is unique.
\begin{thm}\textup{(Main Theorem of Polynomial Interpolation)}\\
Let $(x_1,y_1),(x_2,y_2),...,(x_n,y_n)$ be $n$ points in the plane with distinct $x_i$. There exists one and only one polynomial $P(x)$ of degree less than or equal to $n-1$ satisfying $P(x_i)=y_i \,for \,i = 1,...,n$.
\end{thm}


Because of its special properties, like smoothness, polynomials are often used in interpolation. The Lagrange interpolating polynomial, which has the form 
\begin{align}
P(x)= \sum_{k=1}^n \left( \prod_{j\neq k}\frac{x-x_j}{x_k-x_j} \right) y_k,
\end{align}
is one of the representations of the interpolating polynomial. However, $(2)$ must computed repeatedly each time when adding node points. An alternative way, which has less computational cost, is called Newton interpolation. The Newton interpolating polynomial has the following form:
\begin{align*}
P(x)= &y_1+ f[x_1,x_2](x-x_1)+ f[x_1,x_2,x_3](x-x_1)(x-x_2)\\&+ ... + f[x_1,x_2,...,x_n](x-x_1)(x-x_2)(x-x_3)...(x-x_n-1).
\end{align*}
where $f[x_i,x_j] = \frac{y_i-y_j}{x_i-x_j}$, and $f[x_1,x_2,...,x_k]= \frac{f[x_1,x_2,...,x_k-1]-f[x_2,x_3,...,x_k]}{x_1-x_k}.$  
Both Lagrange interpolation and Newton interpolation require the polynomial to interpolate the given data points, i.e. $P(x_i)= y_i$. Notice that given the same data points $(x_i,y_i)$ for $i = 1,2,...n$, the polynomials derived from Interpolation methods like Newton and Lagrange are the same. Hermite Interpolation requires the polynomial to satisfy the interpolation conditions on derivatives, i.e. to find the polynomial that satisfies 
\begin{align*}
P(x_i) = y_i,\quad P'(x_i) = y'_i,\quad for\, i = 1,2,...,n
\end{align*}
where $y'_i$ are the derivatives of original functions at the node point $(x_i,y_i)$.
\item

Piecewise Interpolation\\
Five interpolation results are shown in figure \ref{F1}. Notice that the original function $f(x)= \frac{1}{x^2+1}$ is shown in black, and the Lagrange interpolating polynomial is shown in red. Figure \ref{F1} shows that a high order interpolating polynomial may have unstable result at the edge of an interval. Piecewise Interpolation is often used to avoid increasing the degree of polynomial when adding more node points. 

Nearest Neighbour Interpolation is the simplest interpolation method in one dimension. The result of this interpolation is a piecewise constant function constructed by selecting the value of nearest given data point and ignoring the neighbour points in each interval. In figure \ref{F1}, the result of nearest neighbour interpolation for approximating $f(x)=\frac{1}{x^2+1},x\in [-5,5]$ are shown in pink. Compared to the original function, this method does not guarantee good accuracy, but it is a fast interpolation method on two dimensional problems. Another simple piecewise interpolation is based on Linear Interpolation. Suppose we have two node points, $(x_1,y_1),(x_2,y_2)$. Then the linear interpolating function is of the form 
\begin{align*}
P(x) = y_1 + (y_2 - y_1) \frac{x -  x_1}{x_2 - x_1} \mbox{ at the point } (x,y).
\end{align*}
In Piecewise Linear Interpolation, one finds a piecewise function $P(x)$ satisfying: 
\begin{enumerate}
\item $P(x) \in C[a,b];$
\item $P(x_i) = y_i(i=1,2,...,n);$
\item $P(x)$ is a linear function in each subinterval $[x_i,x_i+1]$.
\end{enumerate} 

The result of Piecewise Linear Interpolation is shown in the figure \ref{F1}. Apparently, it is not the most accurate method especially when the number of node points is not large enough. Piecewise Cubic Hermite interpolation is shown in the figure \ref{F1} in blue. The main idea of this method, based on Hermite  interpolation, is to construct a continuously differentiable piecewise polynomial $P(x)$. Although Piecewise Linear Interpolation and Piecewise Hermite Interpolation can avoid the problem of oscillation, the functions constructed by using these two methods usually are not smooth curves. The curve constructed by Cubic Spline Interpolation is connected by piecewise cubic curves which are twice continuously differentiable at the connection points.\\

\begin{figure}[H]
\centering
\begin{minipage}[c]{0.8\textwidth}
\centering
\includegraphics[width=10cm]{interpolationdemo}
\caption{Four One Dimension Interpolation Methods}
\label{F1}
\end{minipage}
\end{figure}

\end{enumerate}

\subsection{Multivariate Interpolation} 

If the interpolated function has more than one variable, Multivariate Interpolation such as, Nearest Neighbour Interpolation, Bilinear Interpolation, and Bicubic Interpolation are usually used to approximate it. In Figure \ref{F2}, these three interpolation results for given data sets are drawn on a uniform two dimensional grid. The interpolated function is $f(x,y)= \sin y\cos x$. As mentioned before, Nearest Neighbour Interpolation is a fast interpolation method that can be used in any dimension. Shown in Figure \ref{F2}, the interpolation function is also an piecewise constant function. Bilinear Interpolation is the extension of Linear Interpolation in two dimensions. The main process of Bilinear Interpolation is to do the interpolation first in the $x$ -direction(or $y$ -direction) and then find the estimation by doing interpolating in the $y$ -direction(or $x$ -direction). Bicubic Interpolation can construct a smooth interpolation curve, see also in figure \ref{F2}.


Previous interpolation methods are commonly used in different areas, however, these methods ignored the physical meaning of the given data sets. Also, in real life problems, like rainfall or temperature estimation, the data sets are usually not distributed well. In this paper, we want to use harmonic function to interpolate the given data points. We discussed two situations: first, the data points are given on the boundary of the domain; second, we only know data at several points inside the domain. 
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.9\textwidth}
\centering
\includegraphics[width=10cm]{interpolationdemo2}
\caption{Three Two Dimension Interpolation Methods}
\label{F2}
\end{minipage}
\end{figure}



\section{Properties of Harmonic Function}
In the previous section, we introduced several common interpolation methods. As we mentioned before, some interpolation methods have not considered the physical meaning of given data or its distribution. Here, we want to use harmonic functions, one of the most beautiful type of functions to interpolate the given data points. First, let's define what harmonic function is.\\
\begin{Def}
A harmonic function is a twice continuously differentiable function $u:\Omega\to \mathbb{R}$ that satisfies Laplace's equation  where $\Omega$ is the domain in $\mathbb{R}^n$. That is 
\begin{align}
\Delta u = 0
\end{align}
\label{defHarmonic}
\end{Def}
In other words, a harmonic function is a solution of Laplace equation. Some properties of harmonic function will be introduced in the following. It is these properties which make it possible to be considered as one kind of interpolation function. 
In this paper, we fix $\Omega \subset \mathbb{R}^2$ be a bounded smooth domain (open, connected set). Recall
\begin{thm}\textup{(Divergence Theorem)}
Let $\Omega \subset\mathbb{R}^2$, and $F \in C^1(\overline \Omega)$, then we have:
\[
\int_{\Omega}\nabla\cdot F\,dx = \int_{\partial\Omega}F\cdot n\,dS
\]
where $n$ is the unit outward normal.
\label{DT}
\end{thm}

From the above theorem, one can easily derive the following 


\begin{thm}\textup{(Green's First Identity)}
Let $\Omega\subset\mathbb{R}^2$ and  $u \in C^2(\overline \Omega)$, and let $v \in C^1(\overline \Omega)$. Then we have 
\begin{align}
\int_{\Omega}\nabla v\cdot \nabla u \,dx + \int_{\Omega} v\,\Delta u \,dx = \int_{\partial\Omega}v \frac{\partial u}{\partial n}\,dS.
\end{align}
\label{GFI}
\end{thm}
\begin{proof}
From the equality, 
\[
\nabla v\cdot\nabla u + v\,\Delta u = \nabla\cdot(v\nabla u),
\]
and Theorem \ref{DT},we have
\[
\int_{\Omega}\nabla v\cdot\nabla u\,dx+\int_\Omega v\Delta u\,dx
=\int_{\Omega}\nabla\cdot(v\nabla u)dx
=\int_{\partial\Omega}(v\nabla u)\cdot n\,dS
=\int_{\partial\Omega}v\,\frac{\partial u}{\partial n}\,dS.
\]
By symmetry, (4) can be rewritten as the following:
\begin{align}
\int_{\Omega}\nabla v\cdot \nabla u \,dx + \int_{\Omega} u\,\Delta v \,dx = \int_{\partial\Omega}u \frac{\partial v}{\partial n}\,dS.
\end{align}
\end{proof}



\noindent  Let  $g$ be a function  defined on $\partial\Omega$. Define
\[ 
\mathcal{A} = \{ w \in C^2(\Omega)\,|\, w \equiv g  \textup{ on } \partial\Omega\}
\]

\begin{Def}
For $w \in \mathcal {A}$, define the energy functional of $w$ to be 
\begin{align}
I[w]=\frac{1}{2}\int_\Omega \vert \nabla w \rvert^2.
\end{align}
\label{defenergy}
\end{Def}

\begin{thm}\textup{(Dirichlet's Principle)}

\noindent For $u \in \mathcal{A}$, we have
\[
u \textup{ is a harmonic function} \textup{ if and only if } I[u] = \min_{w \in \mathcal{A}} I[w].
\]
\label{DP}
\end{thm}
\begin{proof}
We first prove necessity. For fixed $v \in C_{c}^\infty(\Omega)$, we can define $\phi(t) = I[u + tv]$, note that $u+tv \in A$ because $v$ has compact support. Note that 
\begin{align*}
f(t)
&= \frac{1}{2}\int_\Omega |\nabla (u+tv)|^2\,dx \\ 
&= \frac{1}{2} \int_\Omega |\nabla u + t\nabla v|^2\,dx \\
&= \frac{1}{2} \int_\Omega [\,|\nabla u |^2+ 2t\,\nabla u\nabla v + t^2\, |\nabla v|^2]\,dx \\
&= I[u] + t\,\int_\Omega \nabla u \nabla v + t^2\, I[v].
\end{align*}
Based on our assumption, $f(u)$ is minimal at $t=0$. Hence, consider its derivative with respect to t
\begin{align*}
0&= \frac{d}{dt}f(t)|_{t=0}= \int_\Omega \nabla u \nabla v \,dx \\
&= \int_{\partial\Omega} v\, \frac{\partial u}{\partial n}\,dS - \int_{\Omega} v\, \Delta u\, dx.
\end{align*}
Since $v$ has compact support, therefore $v = 0$ on $\partial\Omega$, and it follows that
\[
\int_{\Omega} v\, \Delta u\, dx = \, 0    \,\, \,\,\forall v \in C^2_{c}(\Omega).
\]

Here we claim that $\Delta u =\, 0$. Suppose $\Delta u(x) \neq 0$ for some $x \in \Omega$, we may assume $\Delta u(x) > 0$, also notice that $\Delta u$ is continuous at $x$, therefore, there exists $B(x,r) \subset \Omega$ such that $u(y)>0, \,\,\textup{ for all } y \in B(x,r)$. We can also pick $v \in C^2_{c}(B(x,r))$ such that $v(y) > 0\,\,\,for\,\, y\in B(x,r)$.\,Then, we have:\\
\[
\int_{\Omega} v\Delta u \,dx = \int_{B(x,r)} v\Delta u \,dx \,>\,0.
\]
Contradiction. 
\par

To prove the other direction, let $w \in A$, and we need to show that $I[u] \le I[w]$. \\
Let $v = w - u$, then we have: \\
\begin{align*}
I[w]
&= I[v + u] = \int_{\Omega}|\nabla(v+u)|^2\\
&= \int_{\Omega}(\nabla v + \nabla u)\cdot(\nabla v + \nabla u)\\
&= \int_{\Omega}|\nabla v|^2 +\int_{\Omega}|\nabla u|^2 + 2\int_{\Omega}(\nabla v\,\cdot\nabla u)\\
&= I[v] + I[u] +  \int_{\Omega}\nabla(w-u)\cdot\nabla u \\
&= I[v]\, + I[u] + \, \int_{\Omega}(w-u)\cdot\Delta u -\int_{\partial\Omega}(w-u)\frac{\partial u}{\partial n}dS\\
& = I[v] + I[u] \ge I[u].
\end{align*}
\end{proof}
The Dirichlet Principle states that if $u$ is a harmonic function that satisfies $(3)$, then it minimizes the energy.
To find the approximation function for the given data or image, we usually want to find the minimizer or at least the local minimizer of some functional. By the Dirichlet Principle, the harmonic function is the minimizer of the gradient $u$ in $L^2$ norm. Before using harmonic functions to interpolate the given Dirichlet boundary data, one question is why only one solution exists. The uniqueness can be explained by the Maximum Principle. I will prove it by using the following lemma called the Mean Value Property. 

\begin{lemma}\textup{(Mean Value Property)}
Suppose $u$ is a harmonic function, and $u \in C^2(\Omega)$. For $x \in \Omega$ and any open ball $B(x,r) \subset \Omega$, we have
\begin{align*}
u(x)
&= \fint_{B(x,r)} u\,dy =\frac{1}{|B(x,r)|} \int_{B(x,r)} u\,dy\\
&= \fint_{\partial B(x,r)} u\,dy = \frac{1}{|\partial B(x,r)|} \int_{\partial B(x,r)} u\,dy.
\end{align*}
\label{MVP}
\end{lemma}

\begin{thm}\textup{(Maximum Principle)}\\
If $u$ is a harmonic function in $\Omega$, and $u \in C^2(\Omega)\cup C(\partial\Omega)$, then\\
\begin{enumerate}
\item we have \[
\max \limits_{\overline{\Omega}} u = \max \limits_{\partial\Omega}u.
\]
\item If there exist $ x_{0} \in \Omega$ such that $u(x_{0}) = \max \limits_{\partial\Omega}u$, and $\Omega$ is connected, then $u \equiv u(x_{0})$ in $\Omega$.
\end{enumerate}
\label{MP}
\end{thm}
Suppose we have two harmonic functions $u_1,u_2$, and $u_1,u_2 \in C^2(\Omega)\cup C(\partial\Omega)$. Suppose both of them satisfy the given boundary data. Consider the function $ u_1 - u_2$, it has to be zero on the boundary $\Omega$. By the first part of the Maximal Principal, we know that $u_1 - u_2 \leq 0$ on $\Omega \ \partial \Omega$. Now we switch order, and consider function $u_2 - u_1$. Apply Maximum Principal again we get $u_2 - u_1 \leq 0$. Hence, we know that $u_1 = u_2$ on $\Omega$. 

\begin{thm}\textup{(Green's Second Identity)}
Suppose $u,v\in C^2(\Omega)\cup C^1(\overline{\Omega})$. Then,
\[
\int_{\Omega}u\Delta v -v\Delta u\,dx= \int_{\partial\Omega}u\frac{\partial v}{\partial n} - v\frac{\partial u}{\partial n}\,dS.
\]
Here, $\partial/\partial n$ means the directional derivative with respect to the unit outer normal.
\label{GSI}
\end{thm}

\begin{proof}
In view of Green's First Identity, we have
\[
\int_{\Omega}\nabla u\nabla v\,dx + \int_{\Omega} u\Delta v\,dx = \int_{\partial\Omega} u\frac{\partial v}{\partial n}\,dS
\]
and
\[
\int_{\Omega}\nabla v\nabla u\,dx + \int_{\Omega} v\Delta u\,dx = \int_{\partial\Omega} v\frac{\partial u}{\partial n}\,dS.
\]
The result follows from subtracting these two equations.
\end{proof}

In the reminder of this thesis, we focus our attention on $n=2$.

\begin{Def}
Define
\[
\Phi(x)=-\frac{1}{2\pi}\log|x|
\]
for $x\in\mathbb{R}^2\backslash\{0\}$. The function $\Phi$ is called \textbf{the fundamental solution of Laplace equation} (in dimension $2$). We also define $\Phi_x(y) = \Phi(x-y)$.
\label{FSLE}
\end{Def}

It is easy to see that $\Delta_y \Phi_x(y)=0$ for $y \ne x$ . Now we can state our main theorem.

\begin{thm}
Let $u\in C^2(\Omega)\cup C^1(\overline{\Omega})$ be harmonic.
\begin{enumerate}
\item For $x\in\Omega$, we have
\[
u(x) = \int_{\partial\Omega}\Phi(y-x)\frac{\partial u}{\partial n}(y)-u(y)\frac{\partial\Phi}{\partial n}(y-x)\,dS(y)
\]
\item Suppose $\partial\Omega$ is smooth. For $x\in \partial\Omega$, we have
\[
\frac{1}{2}u(x)
=
\int_{\partial\Omega}\Phi(y-x)\frac{\partial u}{\partial n}(y)-u(y)\frac{\partial\Phi}{\partial n}(y-x)\,dS(y);
\]
\end{enumerate}
\label{BoundaryIntergral}
\end{thm}

\begin{proof}
Suppose $x\in\Omega$. Since $\Omega$ is open, we can find $\epsilon>0$ such that $\overline{B_{2\epsilon}(x)}\subset\Omega$. Let $\Omega_\epsilon = \Omega\backslash B_\epsilon(x)$ Now we apply Green's Second Identity to obtain
\[
0=\int_{\partial\Omega_\epsilon}u(y)\frac{\partial\Phi}{\partial n}(y-x)-\Phi(y-x)\frac{\partial u}{\partial n}(y)\,dS(y).
\]
Because $\partial\Omega_\epsilon=\partial\Omega\cup\partial B_\epsilon(x)$, we have
\[
\int_{\partial B_\epsilon(x)}u(y)\frac{\partial\Phi}{\partial n}(y-x)-\Phi(y-x)\frac{\partial u}{\partial n}(y)\,dS(y)
=
\int_{\partial\Omega}\Phi(y-x)\frac{\partial u}{\partial n}(y)-u(y)\frac{\partial\Phi}{\partial n}(y-x)\,dS(y).
\]
Our aim is to show that the left-hand-side will converge to $u(x)$ as $\epsilon\to 0$. From the Mean Value Property, we obtain
\begin{align*}
\int_{\partial B_\epsilon(x)}u(y)\frac{\partial\Phi}{\partial n}(y-x)\,dS(y)
&=\int_{\partial B_\epsilon(x)}u(y)\frac{1}{2\pi\epsilon}\,dS(y)=u(x).
\end{align*}
Remember that the harmonic function is smooth and that $\overline{B_{2\epsilon}(x)}$ is compact. For each $y\in\overline{B_\epsilon(x)}$ we have
\eat{\[
\left|\frac{\partial u}{\partial n}(y)\right|
\leq
\frac{2^{2+1}\cdot 2}{\pi\epsilon^{2+1}}||u||_{L^1(B_\epsilon(y))}
\leq
\frac{16}{\pi\epsilon^3}||u||_{L^1(\overline{B_{2\epsilon}(x)})}
\triangleq M<\infty.
\]}
\[
\left|\frac{\partial u}{\partial n}(y)\right|
\leq
\left|\nabla u(y)\right|
\leq
M,
\]
where $M$ is the maximum value of $|\nabla u|$ on $\overline{B_{2\epsilon}(x)}$. Hence, we get
\begin{align*}
\left|\int_{\partial B_\epsilon(x)}\Phi(y-x)\frac{\partial u}{\partial n}(y)\,dS(y)\right|
&\leq
\int_{\partial B_\epsilon(x)}\left|\Phi(y-x)\frac{\partial u}{\partial n}(y)\right|\,dS(y) \\
&\leq
M\int_{\partial B_\epsilon(x)}|\Phi(y-x)|dS(y)\\
&=
M\cdot\frac{\log\epsilon}{2\pi}\cdot 2\pi\epsilon
=M(\epsilon\log\epsilon).
\end{align*} 

Because $M$ is a uniform bound for all $\epsilon'\leq\epsilon$ and because $\epsilon\log\epsilon\to 0$ as $\epsilon \to 0$, this integral approaches $0$ as $\epsilon\to 0$. Hence, the first part of the theorem follows.

Now suppose $\partial\Omega$ is smooth and let $x\in\partial\Omega$. Pick $\epsilon>0$ such that $B_{2\epsilon}(x)\cap \Omega\subset \Omega$. By similar argument, we derive
\[
\int_{\partial B_\epsilon(x)}u(y)\frac{\partial\Phi}{\partial n}(y-x)-\Phi(y-x)\frac{\partial u}{\partial n}(y)\,dS(y)
=
\int_{\Gamma_\epsilon}\Phi(y-x)\frac{\partial u}{\partial n}(y)-u(y)\frac{\partial\Phi}{\partial n}(y-x)\,dS(y),
\]
where $\Gamma_\epsilon = \partial\Omega_\epsilon\backslash B_\epsilon(x)$. Since $\partial\Omega$ is smooth, the limit of left-hand-side as $\epsilon\to 0$ will be
\[
\lim_{\epsilon\to 0}\int_0^\pi \left[u(\epsilon,\theta)\frac{\partial\Phi}{\partial n}(\epsilon)-\Phi(\epsilon)\frac{\partial u}{\partial n}(\epsilon,\theta)\right]\epsilon\,d\theta.
\]
Let's estimate this integral. The first part is
\[
\int_0^\pi u(\epsilon,\theta)\frac{\partial\Phi}{\partial n}(\epsilon)\epsilon\,d\theta
=\int_0^\pi u(\epsilon,\theta)\frac{1}{2\pi\epsilon}\epsilon\,d\theta
=\frac{1}{2\pi}\int_0^\pi u(\epsilon,\theta)d\theta
\to \frac{1}{2\pi}\pi u(x)
\]
as $\epsilon\to 0$. Using a similar trick to find a uniform bound $M$, the second part is
\[
\left|\int_0^\pi\Phi(\epsilon)\frac{\partial u}{\partial n}(\epsilon,\theta)\epsilon\,d\theta\right|
\leq M\cdot\frac{\log\epsilon}{2\pi}\cdot\pi\epsilon
=\frac{M}{2}\epsilon\log\epsilon
\to 0
\]
as $\epsilon\to 0$. Because $\Gamma_\epsilon\to\partial\Omega$ as $\epsilon\to 0$, we conclude that
\[
\frac{1}{2}u(x)
=
\int_{\partial\Omega}\Phi(y-x)\frac{\partial u}{\partial n}(y)-u(y)\frac{\partial\Phi}{\partial n}(y-x)\,dS(y);
\]
the second part of the theorem follows.
\label{BEM}
\end{proof}




\section{Method and Numerical Experiment}
Our main goal is to find the solution of the following problem:
\begin{align}
u(x_i) = v_i, \, \Delta u = 0,
\end{align}
where $v_i \in \mathbb{R}$ are given function values of an unknown smooth harmonic function which has fixed boundary or interior values at given data points $x_i \in \mathbb{R}^2$. Two different situations are considered here: given interpolating nodes are on the boundary, and when the given interpolating nodes are distributed on the whole domain. 


\subsection{Data are given on Boundary}
\subsubsection{Method}
\eat{Based on the theorem 3.3 we proved in the previous section,  for any given boundary data, we consider function $u$ satisfies $\Delta u = 0$ and $u(x,y) = f(x,y)$ on its boundary.}
Now from both the Dirichlet Principle and Maximum Principle, we consider finding the unique interpolating harmonic function on the given boundary. However, it is not easy to find the harmonic function. Let's consider the following equation system:
\begin{equation}  
\left\{  
             \begin{array}{lr}
             \frac{\partial u}{\partial t} = \Delta u \quad \textup{ for } (x,y,t) \in \Omega \times (0, \infty),\\
             u(x,y,0) = F(x,y) \quad \textup{ for } (x,y) \in \Omega\\
             u(x,y,t) = f(x,y)  \quad \textup{ for } (x,y,t) \in \partial \Omega \times [0, \infty)      
             \end{array}  
\right.
\label{heat}
\end{equation}  

Hence, instead of finding the harmonic function directly, we use the heat equation to find it. Why this is correct? For a function of two variable and the time variable $t$, like the heat equation in (\ref{heat}), we have the following:
\begin{align*}
\lim_{t \to \infty} u(x,y,t) = v(x,y).
\end{align*}
Hence  $v(x,y)$ satisfies:

\begin{align*}  
\left\{  
             \begin{array}{lr}
             \Delta v(x,y) = 0 \quad \textup{ for } (x,y) \in \Omega\\
             v(x,y) = f(x,y)  \quad \textup{ for } (x,y) \in \partial \Omega.     
             \end{array}  
\right.
\end{align*} 
%This is because as $t \to \infty$, the heat conduction goes to steady-state, that is $\frac{\partial u}{\partial t} = 0$. 

To find the approximation by heat equation, a known method called finite differences, which replaces each derivative by a difference quotient, will be used in the experiment.To find the function satisfying (\ref{heat}) on $\Omega$, we do the following steps:
\begin{enumerate}
\item[•]
Draw a grid on $\Omega$ and choose mesh size $\Delta x, \Delta y, \Delta t$ for $x,y,t$. 
\item[•]
Approximate the value of $u(i\Delta x, j\Delta y, n\Delta t)$ for $x = i\Delta x, y = j\Delta y, t = n\Delta t$ by a number $u_{i,j}^n $ indexed by integers i, j, n:

\[
 u_{i,j}^n = u(i\Delta x ,j\Delta y, n\Delta t)  
\] 
Therefore we can immediately have the following approximations:
\[
\frac{\partial u}{\partial t} \thicksim \frac{u_{i,j}^{n+1}-u_{i,j}^{n}}{\Delta t}
\]
\[
\frac{\partial^2 u}{\partial x ^2} \thicksim \frac{u_{i+1,j}^{n}-2u_{i,j}^{n}+ u_{i-1,j}^{n}}{(\Delta x) ^2}
\]
\[
\frac{\partial^2 u}{\partial y ^2} \thicksim \frac{u_{i,j+1}^{n}-2u_{i,j}^{n}+ u_{i,j-1}^{n}}{(\Delta y) ^2}
\]
\item[•]
One can quickly rewrite $\frac{\partial u}{\partial t} = \Delta u$, by using finite differences: 
\begin{align}
\frac{u_{i,j}^{n+1}- u_{i,j}^n}{\Delta t}
 = \frac{u_{i+1,j}^{n}-2u_{i,j}^n + u_{i-1,j}^n}{(\Delta x)^2} +  
   \frac{u_{i,j+1}^{n}-2u_{i,j}^n + u_{i,j-1}^n}{(\Delta y)^2}
 \label{NE11}
\end{align}
Drawing a $(N+1) \times (N+1)$ grid on $\Omega$, and choose the same mesh size for both $x$ and $y$, that is $\Delta x = \Delta y$, (\ref{NE11}) can be simplified to 
\begin{align}
u_{i,j}^{n+1} 
 = \frac{\Delta t}{(\Delta x) ^2}(u_{i+1,j}^{n}+ u_{i-1,j}^n - 4u_{i,j}^{n}+ u_{i,j+1}^{n} + u_{i,j-1}^n) + u_{i,j}^n,
 \label{NE12}
\end{align}
\item[•]
$u_{i,j}^0$ is given by the second equation in (\ref{heat}). Using (\ref{NE12}) gives us $u_{i,j}^1$, then (\ref{NE12}) gives $u_{i,j}^2$ and so on. \\
Finally, (\ref{heat}) can be rewritten numerically as the following:
\begin{equation}  
\left\{  
             \begin{array}{lr}
             u_{i,j}^{n+1} = \frac{\Delta t}{(\Delta x) ^2}(u_{i+1,j}^{n}+ u_{i-1,j}^n - 4u_{i,j}^{n}+ u_{i,j+1}^{n} + u_{i,j-1}^n) + u_{i,j}^n\\
             u_{i,j}^0 = F_{i,j} \quad\textup{ for } i,j = 0,1,2,...N\\
             u_{i,j}^n = f_{i,j} \quad\textup{ for } i,j=0, N,  
             \end{array}  
\right.
\label{FN}
\end{equation}
\end{enumerate} 

 


\subsubsection{Numerical Experiment}
\eat{In this part, we trying to find different types of harmonic function using iteration on time of (\ref{FN}).} In all numerical simulations, we take as the domain $\Omega$ is the square with boundary values given. That is, $\Omega$ is set as $[-0.5,0.5] \times [-0.5,0.5]$ and mesh points are placed on $(N+1) \times (N+1)$  gridding points on it.  As described in equation (\ref{FN}), $u_{i,j}^0$ is given on the whole domain, and the boundary value is fixed. In the following tables, variable $k$ is the number of iterations, $\Delta x$ is the mesh size of chosen domain, $\frac{\Delta t}{(\Delta x)^2}$ is the step size of each iteration. In all numerical simulation, step size is $0.25$.  The errors shown are the investigation of accuracy compared to the original harmonic function on $\Omega$. 

\begin{enumerate}
\item[(1)]
$v(x,y) = x^2 - y^2 + xy$
\begin{table}[htb]
\caption{Experiment Report for Function (1)}
\centering
\setlength{\tabcolsep}{5mm}{
\begin{tabular}{c|c|c|c|c}
\hline
k & deltx & Max & Mean & Var \\
\hline
$2500$ & $0.025$ & $3.594\times10^{-15}$  & $1.383\times10^{-15}$ &$1.162\times10^{-30}$ \\
\hline
$5000$& $0.025$ & $4.718\times10^{-16}$ & $1.629\times10^{-16}$&$1.718\times10^{-32}$ \\
\hline
$2500$ & $0.0125$ & $4.583\times10^{-5}$ & $1.812\times10^{-5}$ & $1.846\times10^{-10}$\\
\hline
$5000$ & $0.0125$ & $2.000\times10^{-8}$ & $7.906\times10^{-9}$&$3.510\times10^{-17}$ \\
\hline
\end{tabular}}
\end{table}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat1-1-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat1-1-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat1-1-3.pdf}
\end{minipage}
\end{figure}
\vspace{-6cm} 
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat1-2-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat1-2-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat1-2-3.pdf}
\end{minipage}
\end{figure}
\vspace{-6cm} 
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat1-3-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat1-3-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat1-3-3.pdf}
\end{minipage}
\end{figure}
\vspace{-6cm}
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat1-4-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat1-4-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat1-4-3.pdf}
\end{minipage}
\end{figure}
\item[(2)]
$v(x,y) = \frac{x+1}{(x+1)^2+(y+1)^2}$
\begin{table}[htb]
\caption{Experiment Report for Function (2)}
\centering
\setlength{\tabcolsep}{5mm}{
\begin{tabular}{c|c|c|c|c}
\hline
k & deltx & Max & Mean & Var \\
\hline
$5000$ & $0.025$ & $3.594\times10^{-5}$  & $1.160\times10^{-5}$ &$1.170\times10^{-10}$ \\
\hline
$10000$& $0.025$ & $2.785\times10^{-5}$ & $1.153\times10^{-5}$&$1.163\times10^{-10}$ \\
\hline
$5000$ & $0.0125$ & $0.072$ & $0.0068$ & $2.580\times10^{-5}$\\
\hline
$10000$ & $0.0125$ & $3.587\times10^{-4}$ & $1.463\times10^{-4}$&$1.191\times10^{-8}$ \\
\hline
\end{tabular}}
\end{table}
\vspace{-2cm}
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat2-1-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat2-1-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat2-1-3.pdf}
\end{minipage}
\end{figure}
\vspace{-6cm} 
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat2-2-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat2-2-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat2-2-3.pdf}
\end{minipage}
\end{figure}
\vspace{-6cm} 
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat2-3-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat2-3-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat2-3-3.pdf}
\end{minipage}
\end{figure}
\vspace{-6cm}
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat2-4-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat2-4-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat2-4-3.pdf}
\end{minipage}
\end{figure}

\item[(3)]
$v(x,y) = \log (\sqrt{(x+1)^2+(y-1)^2})$
\begin{table}[htb]
\caption{Experiment Report for Function (3)}
\centering
\setlength{\tabcolsep}{5mm}{
\begin{tabular}{c|c|c|c|c}
\hline
k & deltx & Max & Mean & Var \\
\hline
$5000$ & $0.025$ & $1.454\times10^{-5}$  & $5.425\times10^{-6}$ &$1.813\times10^{-11}$ \\
\hline
$10000$& $0.025$ & $1.461\times10^{-5}$ & $5.467\times10^{-6}$&$1.836\times10^{-11}$ \\
\hline
$5000$ & $0.0125$ & $0.0119$ & $0.0047$ & $1.234\times10^{-5}$\\
\hline
$10000$ & $0.0125$ & $2.481\times10^{-4}$ & $9.780\times10^{-5}$&$5.396\times10^{-9}$ \\
\hline
\end{tabular}}
\end{table}
\vspace{-2cm}
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat3-1-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat3-1-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat3-1-3.pdf}
\end{minipage}
\end{figure}
\vspace{-6cm} 
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat3-2-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat3-2-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat3-2-3.pdf}
\end{minipage}
\end{figure}
\vspace{-6cm} 
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat3-3-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat3-3-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat3-3-3.pdf}
\end{minipage}
\end{figure}
\vspace{-6cm}
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat3-4-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat3-4-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat3-4-3.pdf}
\end{minipage}
\end{figure}


\item[(4)]
$v(x,y) =e^{x}(x\cos{y}-y\sin{y})$
\begin{table}[H]
\caption{Experiment Report for Function (4)}
\centering
\setlength{\tabcolsep}{5mm}{
\begin{tabular}{c|c|c|c|c}
\hline
k & deltx & Max & Mean & Var \\
\hline
$2000$ & $0.025$ & $3.342\times10^{-5}$  & $1.466\times10^{-5}$ &$1.031\times10^{-10}$
\\
\hline
$5000$& $0.025$ & $3.156\times10^{-5}$ & $1.390\times10^{-5}$&$9.203\times10^{-11}$ \\
\hline
$2000$ & $0.0125$ & $0.0085$ & $0.0033$ & $6.134\times10^{-6}$\\
\hline
$5000$ & $0.0125$ & $4.682\times10^{-5}$ & $1.303\times10^{-5}$&$1.837\times10^{-10}$ \\
\hline
\end{tabular}}
\end{table}
\vspace{-2cm}
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat4-1-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat4-1-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat4-1-3.pdf}
\end{minipage}
\end{figure}
\vspace{-6cm} 
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat4-2-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat4-2-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat4-2-3.pdf}
\end{minipage}
\end{figure}
\vspace{-6cm} 
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat4-3-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat4-3-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat4-3-3.pdf}
\end{minipage}
\end{figure}
\vspace{-6cm}
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat4-4-1.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat4-4-2.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=6cm]{mat4-4-3.pdf}
\end{minipage}
\end{figure}

\end{enumerate}

\subsection{Scattered Data are given on whole domain}
\subsubsection{Method}
Sometimes, data are not only given on boundaries, but can be given on the whole domain, and not distributed well. Here we use the Boundary Element Method to find the Harmonic function. The theorem (\ref{BEM}) proved in the previous section tells us a harmonic function $u \in C^2(\Omega)\cup C^1(\overline{\Omega})$ can be described by the boundary integral equation: 
\begin{align}
\alpha(x) u(x) = \int \Phi_x \frac{\partial u}{\partial n} - u\frac{\partial \Phi_{x}}{\partial n}dS,
\label{BME2}
\end{align}
where $\alpha(x) = 1 \textup{ for } x \in \Omega \textup{ and } \alpha(x) =\frac{1}{2} \textup{ for } x \in \partial\Omega$.

Now let's replace $\partial \Omega$ by some straight line segments $\Gamma_1, \Gamma_2,...\Gamma_n$. Also, let $y_1,y_2,...y_n \in \mathbb{R}^2$ be the middle points of $\Gamma_1, \Gamma_2,...\Gamma_n$. 
Therefore the approximation of (\ref{BME2}) is: 
\begin{align}
\alpha(x) u(x) 
&= \int \Phi_{x} \frac{\partial u}{\partial n} - u\frac{\partial \Phi_{x}}{\partial n}dS \nonumber \\
&\approx \int_{\Gamma_1 \cup \Gamma_2,...\cup \Gamma_n} \Phi_x
 \frac{\partial u}{\partial n} - u\frac{\partial \Phi_{x}}{\partial n}dS \nonumber \\
&= \sum_{j=1}^{n} \int _{\Gamma_j} \Phi_x \frac{\partial u}{\partial n} - \sum_{j=1}^{n} \int_{\Gamma_j} u \frac{\partial \Phi_x}{\partial n} 
\label{Ds1}
\end{align}

If we use the value of $u(y_i)$ to represent the value of $u$ on each $\Gamma_i$, then (\ref{Ds1}) can be approximated as:
\begin{align*}
\sum_{j=1}^{n} \frac{\partial u(y_j)}{\partial n} \int_{\Gamma_j} \Phi_{x_i}dS- \sum_{j=1}^{n} u(y_j) \int_{\Gamma_j} \frac{\partial \Phi_{x}}{\partial n}
\end{align*}

Hence if we know the value of $u(y_j)$ and value of $\frac{\partial u(y_j)}{\partial n}$ for $i=1,2,...n$ then for each $x \in \overline{\Omega}$ the value of $u(x)$ can be approximated.

Now for given data $x_1,x_2,...x_n \in \mathbb{R}^2$, and $v_1,v_2,...v_n \in \mathbb{R}$, to find $u \in C^2(\mathbb{R}^2)$ such that $\Delta u =0, u(x_i) = v_i $ for all $i$, we do the following: 
\begin{enumerate}
\item[•]
Enclose $x_1,x_2,...x_n$ by a rectangle $\Omega$.
\item[•]
Discretize $\partial \Omega$ into straight line segments $\Gamma_1,\Gamma_2,...\Gamma_n$ and define $y_j$ be the midpoints of $\Gamma_j$.
\item[•]
For each $x_i \in R$ we consider:
\begin{align*}
v_i = \sum_{j=1}^{n} \frac{\partial u}{\partial n} (y_j) \int_{\Gamma_j} \Phi_{x_i}dS - \sum_{j=1}^{n} u(y_j)\int_{\Gamma_j} \frac{\partial \Phi_{x_i}}{\partial n}dS,
\end{align*}
and
\begin{align*}
\frac{1}{2}u(y_i) = \sum_{j=1}^{n} \frac{\partial u}{\partial n} (y_j)\int_{\Gamma_j} \Phi_{x_i}dS - \sum_{j=1}^{n} u(y_j)\int_{\Gamma_j} \frac{\partial \Phi_{x_i}}{\partial n}dS,
\end{align*}

In terms of matrix multiplication we have:\\
\begin{align}
\left( \begin{matrix}
	\int_{\Gamma_j}\Phi_{x_i}dS
\end{matrix}\right)
\left( \begin{matrix}
	\frac{\partial u}{\partial n}(y_1) \\
	\vdots\\
	\frac{\partial u}{\partial n}(y_n)
\end{matrix}\right)  - 
\left( \begin{matrix}
	\int_{\Gamma_j} \frac{\partial \Phi_{x_i}}{\partial n}dS
\end{matrix}\right)
\left( \begin{matrix}
	u(y_1) \\
	\vdots\\
	u(y_n)
\end{matrix}\right) =
\left( \begin{matrix}
	v(1)\\
	\vdots\\
	v(n)
\end{matrix}\right)
\label{MX1}
\end{align} 
and,
\begin{align}
\left( \begin{matrix}
	\int_{\Gamma_j}\Phi_{y_i}dS
\end{matrix}\right)
\left( \begin{matrix}
	\frac{\partial u}{\partial n}(y_1) \\
	\vdots\\
	\frac{\partial u}{\partial n}(y_n)
\end{matrix}\right)  - 
\left( \begin{matrix}
	\int_{\Gamma_j} \frac{\partial \Phi_{y_i}}{\partial n}dS + \frac{1}{2}\delta_{i,j}
\end{matrix}\right)
\left( \begin{matrix}
	u(y_1) \\
	\vdots\\
	u(y_n)
\end{matrix}\right) =
0,
\label{MX2}
\end{align}
where 
\begin{equation*}  
\delta_{i,j}=
\left\{  
             \begin{array}{lr}
             1 \quad if \quad i = j\\
             0 \quad if \quad i \neq j                    
             \end{array}  
\right..
\end{equation*}
Let's combine (\ref{MX1}) and (\ref{MX2}), that is
\begin{align}
	\left(
	\begin{array}{c|c}
	\int_{\Gamma_j}\Phi_{x_i}dS & - \int_{\Gamma_j} \frac{\partial \Phi_{x_i}}{\partial n}dS \\ \hline	
	\int_{\Gamma_j}\Phi_{y_i}dS & - \int_{\Gamma_j} \frac{\partial \Phi_{y_i}}{\partial n}dS - \frac{1}{2}\delta_{i,j}
	\end{array}
	\right)
	\left(
	\begin{array}{c}
	\frac{\partial u}{\partial n}(y_1) \\ 
	\vdots \\
	\frac{\partial u}{\partial n}(y_n)\\  \hline
	u(y_1) \\
	\vdots \\
	u(y_n)
	\end{array}
	\right)=
	\left(
	\begin{array}{c}
	v_1 \\ 
	\vdots\\
	v_n\\ \hline
	0\\
	\vdots\\
	0
	\end{array}
	\right)
\label{MX3}
\end{align}

As long as $u(y_1),...u(y_n)$ and $\frac{\partial u}{\partial n}(y_1),\frac{\partial u}{\partial n}(y_2),...\frac{\partial u}{\partial n}(y_n)$ are known, for each $x_i \in R$, $u(x_i)$ can be computed by 
\begin{align}
\alpha_{x_i} u(x_i) = 
\left(
	\begin{array}{c|c}
	\int_{\Gamma_j}\Phi_{x_i} dS & - \int_{\Gamma_j} \frac{\partial \Phi_{x_i}}{\partial n}dS 
	\end{array}
	\right)
	\left(
	\begin{array}{c}
	\frac{\partial u}{\partial n}(y_1) \\ \hline
	u(y_n) \\
	\end{array}
	\right)
\label{MX4}
\end{align}

\item[•]
If the solution of (\ref{MX4}) is not unique, more constraints are needed here. Remember we proved that lowest energy is attained if it is a harmonic function. Let's say $I(u)$ is energy of the steady field. Hence consider
\begin{align*}
I(u) = \int_R |\nabla u|^2 = \int_{\partial R} u\frac{\partial u}{\partial n}dS,
\end{align*}
also,
\begin{align*}
\int_{\partial R}\frac{\partial u}{\partial n} = 0.
\end{align*}
Discretizing these two equations, we get 
\begin{align*}
I(u) = 
\sum_{j=1}^{n} u(y_j)\frac{\partial u}{\partial n}(y_j)l(\Gamma_j)
\end{align*}
and,
\begin{align*}
\sum_{j=1}^{n}\frac{\partial u}{\partial n}(y_j)l(\Gamma_j)=0
\end{align*}

where $l(\Gamma_i)$ is the length of the $i$-th segment. Hence, we try to solve 
\begin{align*}
\min \{ \sum_{j=1}^{n} u(y_j)\frac{\partial u}{\partial n}(y_j)l(\Gamma_j) \}
\end{align*}
with,
\begin{equation}  
\left\{  
             \begin{array}{lr}
             \left(
	\begin{array}{c|c}
	\int_{\Gamma_j}\Phi_{x_i}dS & - \int_{\Gamma_j} \frac{\partial \Phi_{x_i}}{\partial n}dS \\ \hline	
	\int_{\Gamma_j}\Phi_{y_i}dS & - \int_{\Gamma_j} \frac{\partial \Phi_{y_i}}{\partial n}dS
	\end{array}
	\right)
	\left(
	\begin{array}{c}
	\frac{\partial u}{\partial n}(y_1) \\ 
	\vdots \\
	\frac{\partial u}{\partial n}(y_1)\\  \hline
	u(y_1) \\
	\vdots \\
	u(y_n)
	\end{array}
	\right)=
	\left(
	\begin{array}{c}
	v_1 \\ 
	\vdots\\
	v_n\\ \hline
	0\\
	\vdots\\
	0
	\end{array}
	\right)\\ 
\sum\limits_{j=1}^{n} \frac{\partial u}{\partial n}(y_j)l(\Gamma_j)=0.
    \end{array}
\right.
\label{MX5}
\end{equation}
Notice that for (\ref{MX5}), we do not solve for the exact solution, but a solution such that 
\begin{align}
\left|
\left(
	\begin{array}{c|c}
	\int_{\Gamma_j}\Phi_x dS & - \int_{\Gamma_j} \frac{\partial \Phi_x}{\partial n}dS 
	\end{array}
	\right)
	\left(
	\begin{array}{c}
	\frac{\partial u}{\partial n}(y_1) \\ \hline
	u(y_1) \\
	\end{array}
	\right)-
	\left(
	\begin{array}{c}
	v_1\\
	\vdots\\
	v_n
	\end{array}
	\right)
\right| <\epsilon
\end{align}
for chosen $\epsilon$.
\end{enumerate}  

\subsubsection{Numerical Experiment}
In this part we give the data sets from different harmonic functions and use our method to find the accuracy of proposed scheme. For each simulation, data set $(x_1,v_1),...(x_n,v_n)$ is given on the $N \times N$ grid on $[-0.5,0.5] \times [-0.5,0.5]$. We enclose all the given data points in a bigger rectangular $\Omega$. Then the interpolating result will be compared with the original values of the harmonic function on a $100\times 100$ grid on $\Omega$. The following are the result of the investigation of accuracy compared to the original harmonic function $v$ on $\Omega$.

\begin{enumerate}
\item[(1)]
$v(x,y) = x^2 - y^2 + xy$
\begin{table}[htb]
\caption{Experiment Report for Function (1)}
\centering
\setlength{\tabcolsep}{3mm}{
\begin{tabular}{c|cccc}
\hline
N & 4 & 8 & 16 & 32 \\
\hline
Max of Error & $4.328\times 10^{-2}$ & $3.337\times10^{-2}$  & $1.195\times10^{-2}$ &$6.891\times10^{-3}$ \\
\hline
Min of Error & $5.126\times10^{-7}$ & $6.251\times10^{-10}$ & $1.055\times10^{-12}$&$5.557\times10^{-14}$ \\
\hline
Mean of Error & $4.514\times10^{-3}$ & $1.062\times10^{-3}$ & $1.110\times10^{-4}$ & $2.568\times10^{-5}$\\
\hline
Variance of Error & $3.803\times10^{-5}$ & $8.667\times10^{-6}$ & $4.561\times10^{-7}$&$7.194\times10^{-8}$ \\
\hline
\end{tabular}}
\end{table}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{fig_u_4.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{fig_error_4.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{inter_N=4.pdf}
\end{minipage}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{fig_u_8.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{fig_error_8.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{inter_N=8.pdf}
\end{minipage}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{fig_u_16.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{fig_error_16.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{inter_N=16.pdf}
\end{minipage}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{fig_u_16.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{fig_error_16.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{inter_N=16.pdf}
\end{minipage}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{fig_u_32.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{fig_error_32.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{inter_N=32.pdf}
\end{minipage}
\end{figure}
\item[(2)]
$v(x,y) = \frac{x+1}{(x+1)^2+(y+1)^2}$

\begin{table}[H]
\caption{Experiment Report for Function (2)}
\centering
\setlength{\tabcolsep}{3mm}{
\begin{tabular}{c|cccc}
\hline
N & 4 & 8 & 16 & 32 \\
\hline
Max of Error & $3.294\times 10^{-2}$ & $3.453\times10^{-2}$  & $1.224\times10^{-2}$ & $7.619\times10^{-3}$\\
\hline
Min of Error & $2.491\times10^{-7}$ & $1.852\times10^{-10}$ & $6.239\times10^{-16}$ & $8.882\times10^{-16}$ \\
\hline
Mean of Error & $2.574\times10^{-3}$ & $6.453\times10^{-4}$ & $3.918\times10^{-5}$& $9.288\times10^{-6}$ \\
\hline
Variance of Error & $1.851\times10^{-5}$ & $5.781\times10^{-6}$ & $1.280\times10^{-8}$ & $2.010\times10^{-8}$\\
\hline
\end{tabular}}
\end{table}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{2fig_u_4.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{2fig_error_4.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{2inter_N=4.pdf}
\end{minipage}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{2fig_u_8.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{2fig_error_8.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{2inter_N=8.pdf}
\end{minipage}
\end{figure}


\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{2fig_u_16.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{2fig_error_16.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{2inter_N=16.pdf}
\end{minipage}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{2fig_u_32.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{2fig_error_32.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{2inter_N=32.pdf}
\end{minipage}
\end{figure}


\item[(3)]
$v(x,y) = ln[(x-1)^2+(y-1)^2]$

\begin{table}[H]
\caption{Experiment Report for Function (3)}
\centering
\setlength{\tabcolsep}{3mm}{
\begin{tabular}{c|cccc}
\hline
N & 4 & 8 & 16 & 32 \\
\hline
Max of Error & $9.915\times 10^{-2}$ & $3.406\times10^{-2}$  & $9.309\times10^{-3}$ & $1.1219\times10^{-2}$\\
\hline
Min of Error & $8.743\times10^{-7}$ & $1.004\times10^{-11}$ & $2.487\times10^{-14}$ & $5.863\times10^{-13}$ \\
\hline
Mean of Error & $7.969\times10^{-3}$ & $5.699\times10^{-4}$ & $4.437\times10^{-5}$ & $2.553\times10^{-5}$\\
\hline
Variance of Error & $1.568\times10^{-4}$ & $4.134\times10^{-6}$ & $1.423\times10^{-7}$ & $1.1216\times10^{-7}$\\
\hline
\end{tabular}}
\end{table}
\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{3fig_u_4.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{3fig_error_4.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{3inter_N=4.pdf}
\end{minipage}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{3fig_u_8.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{3fig_error_8.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{3inter_N=8.pdf}
\end{minipage}
\end{figure}


\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{3fig_u_16.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{3fig_error_16.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{3inter_N=16.pdf}
\end{minipage}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{3fig_u_32.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{3fig_error_32.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{3inter_N=32.pdf}
\end{minipage}
\end{figure}


\item[(4)]
$v(x,y) =e^{x}(x\cos{y}-y\sin{y})$
\begin{table}[H]
\caption{Experiment Report for Function (4)}
\centering
\setlength{\tabcolsep}{3mm}{
\begin{tabular}{c|cccc}
\hline
N & 4 & 8 & 16 & 32 \\
\hline
Max of Error & $2.031\times10^{-2}$ & $7.366\times10^{-2}$ & $2.140\times10^{-2}$ & $1.1219\times10^{-2}$\\
\hline
Min of Error & $6.135\times10^{-8}$ & $3.293\times10^{-9}$& $1.727\times10^{-12}$ & $5.863\times10^{-13}$\\
\hline
Mean of Error & $1.577\times10^{-3}$ & $1.849\times10^{-3}$ & $1.111\times10^{-4}$ & $2.442\times10^{-5}$ \\
\hline
Variance of Error & $5.436\times10^{-6}$ & $4.014\times10^{-5}$ & $7.792\times10^{-7}$ & $1.216\times10^{-7}$ \\
\hline 
\end{tabular}}
\end{table}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{4fig_u_4.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{4fig_error_4.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{4inter_N=4.pdf}
\end{minipage}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{4fig_u_8.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{4fig_error_8.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{4inter_N=8.pdf}
\end{minipage}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{4fig_u_16.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{4fig_error_16.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{4inter_N=16.pdf}
\end{minipage}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{4fig_u_32.pdf}
\end{minipage}
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{4fig_error_32.pdf}
\end{minipage}
\centering
\begin{minipage}[t]{0.3\textwidth}
\centering
\includegraphics[width=5cm]{4inter_N=32.pdf}
\end{minipage}
\end{figure}
\end{enumerate}


\section{Conclusion}
This paper focuses on the interpolation method for given data points on $\mathbb{R}^2$. We state that traditional interpolation methods usually do not consider the physical meaning of the interpolating function. Also consider real life problems like weather forecasting, in which data may not be distributed well, we found that using harmonic function as an interpolating function is a reasonable strategy. Two different interpolation methods are discussed here. The first method deals with the situation where data are given on the boundary. It uses the heat equation to approach the steady state on the domain. The experimental result indicated that after enough iterations, the error between the approximated result and the original function could be relatively small. In this method, data are given on the boundary. The second method is a convolution based method, called the boundary element method. In this situation, data are give on the whole domain but not distributed well. The second experiment shows that the approximated data points fit the original harmonic function well even when fewer data points are given. 
Some traditional interpolation methods in one and two dimensions are summarized in the previous section. There are two main reasons why we consider harmonic function as interpolation functions: 
\begin{enumerate}
\item
Some interpolation methods do not consider the physical meaning of data points.
\item
Some interpolation methods cannot be used for scattered data points or data points given only on the boundary.
\end{enumerate}

Hence, a smooth function, harmonic function defined on a open connected domain $\Omega$ are considered here. 
Some properties of the harmonic function are introduced in the third section. Theorem \ref{GFI} can be showed easily based on Divergence Theorem. Then, using Theorem \ref{GFI}, we proved an important property of harmonic functions. Theorem \ref{DP} states that the harmonic function on a bounded open connected domain is the minimizer of the Dirichlet functional. The uniqueness theorem can be proved directly by using Theorem \ref{MP}. These properties are the motivations and reasons for using harmonic function to interpolate the given data. To use the boundary element method, we introduced the fundamental solution of Laplace equation and Theorem \ref{BoundaryIntergral}. Theorem \ref{BoundaryIntergral} tells us as that long as $u \in C^2 (\Omega) \cup C^1(\overline{\Omega})$, then the values of $u$ inside the domain on the boundary of domain can be given by a line integral on the boundary related to its fundamental solution. Two numerical methods are introduced step by step in the fourth section. Experiments are done by Matlab and Python. Both numerical experiment try to investigate the accuracy between the approximated values and the values generated by original harmonic functions. The results demonstrated in the fourth section reveal that both methods works if only boundary data or scatted data are given.

In the future, more work will be done on the application of harmonic function interpolation. There are some different types of interpolation used in different subjects. For example, the method could be used in weather forecasting. In recent years, more interdisciplinary methods have been applied to the study of interpolation methods. In \citep{wong_rainfall_2003}, a neural network methods such as ANN(artificial neural network) was used in interpolation models to predict rainfall. In our future work, one of the topics should be how to create the interpolation model combining our method with different statistical techniques in different areas. 

\nocite{*}
\bibliographystyle{plainnat}
\bibliography{harmonicfunction}

\end{document}
