text:
\section{Enhanced error estimate in $L^2$ \label{sec:duality} }

In this section, we provide an estimate of the $L^2$-norm of the discretization error $u-u_h$, by adapting to the present setting the classical Aubin-Nitsche duality argument. To this end, we need two technical results. The first one is a measure of the non-symmetry of the bilinear form $a_h$.

\begin{lemma}\label{lemma:non-symmetry}
It holds $\forall v,w \in V(\tO;\tT)$
\begin{equation}\label{eq:non-symmetry}
a_h(w,v)-a_h(v,w) = (\partial_\tns w, \nabla v \cdot \bs{d})_{0,\tG} - (\partial_\tns v, \nabla w \cdot \bs{d})_{0,\tG} \,.
\end{equation}
\end{lemma}
\begin{proof}
By \eqref{eq:SBM-a-form}  one has
\begin{equation*}
\begin{split}
&a_h(w,v) = (\nabla w, \nabla v)_{0,\tO} - (\partial_\tns w, \tS v)_{0,\tG} + (\partial_\tns w, \nabla v \cdot \bs{d})_{0,\tG}  \\
& \hskip 3.8cm  -  ( \tS w, \partial_\tns v)_{0,\tG} + \gamma \, (h^{-1} \tS w, \tS v)_{0,\tG} \,.
\end{split}
\end{equation*}
Subtracting from this expression the analogous one in which $w$ and $v$ are exchanged, we obtain the result.
\end{proof}

The second technical result replaces the Galerkin orthogonality property of standard interior discretizations.

\begin{lemma}\label{lemma:non-Gal-orth}
It holds
\begin{equation}\label{eq:non-Gal-Orth}
a_h(u-u_h, v_h) = (R_h u, \partial_\tns v_h)_{0,\tG} - \gamma (h^{-1}R_h u, \tS v_h)_{0,\tG}\,, \qquad \forall v_h \in V_h \,.
\end{equation}
\end{lemma}
\begin{proof}
Consider \eqref{eq:nitsche1} and replace therein the expression of $u-\bar{g}$ given by \eqref{eq:u-g}, obtaining
\begin{equation*}
\begin{split}
&(\nabla u, \nabla v_h)_{0,\tO} \!\! - (\partial_\tns u, v_h)_{0,\tG} \!\! -  ( \tS u + R_h u, \partial_\tns v_h)_{0,\tG} \!\! +  \gamma \, (h^{-1} \tS u +R_h u, \tS v_h)_{0,\tG} \\
& \hskip 3cm = (f,v_h)_{0,\tO} -  (\bar{g}, \partial_\tns v_h)_{0,\tG} + \gamma \, (h^{-1} \bar{g}, \tS v_h)_{0,\tG}\,, \quad \forall v_h \in V_h \,.
\end{split}
\end{equation*}
Subtracting from this identity the definition \eqref{eq:SBM-formulation-1} of SBM solution, we obtain the result.
\end{proof}

We wish to estimate
\begin{equation}\label{eq:Riesz}
\Vert u-u_h \Vert_{0,\tO} = \sup_{\tilde{z} \in L^2(\tO)} \frac{(u-u_h,\tilde{z})_{0,\tO}}{\Vert \tilde{z} \Vert_{0,\tO}}\,.
\end{equation}
Given $\tilde{z} \in L^2(\tO)$, let $z \in L^2(\Omega)$ be its extension by $0$ outside $\tO$, and let $w$ be the solution of the auxiliary boundary-value problem
\begin{equation}\label{eq:aux-Dirichlet}
\begin{split}
-\Delta w &= z \quad \text{in \ } \Omega \,, \\
w &=0 \quad \text{on \ } \Gamma \,.
\end{split}
\end{equation}
Let $s_0 \in (\tfrac12,1]$ be such that $w \in H^{1+s_0}(\Omega)$; thus, $s_0=1$ if there are no boundary singularities ($\lambda_S=\infty$ in \eqref{eq:min-lambda}), otherwise it is any real number satisfying $s_0 < \lambda_S<1$. It holds $\Vert w \Vert_{1,\Omega} + |\nabla w|_{s_0,2,\Omega} \leq c \Vert z \Vert_{0, \Omega}= c \Vert \tilde{z} \Vert_{0, \tO}$ for some $c>0$ independent of $\tilde{z}$. Thus, $w \in W(\tO;\tT) := V(\tO;\tT,s_0,2)$ as defined in \eqref{eq:def-Wh1}.

The same argument that led to \eqref{eq:trace-u} shows that the trace $\tilde{q}$ of $w$ on $\tG$ satisfies
\begin{equation}\label{eq:trace-w}
\tilde{q} = -\nabla w \cdot \bs{d} - R(w,\bs{d})\,,
\end{equation}
whence on $\tG$ it holds
\begin{equation}\label{eq:w-q}
0= w-\tilde{q}  \ = \ \tS w  + R_h w \,.
\end{equation}
Thus, for any $v \in V(\tO;\tT)$ we have
\begin{equation*}
\begin{split}
(\tilde{z},v)_{0,\tO} &= -(\Delta w,v )_{0,\tO}  = (\nabla w, \nabla v)_{0,\tO} - (\partial_\tns w, v)_{0,\tG} \\
&= (\nabla w, \nabla v)_{0,\tO} - (\partial_\tns w, v)_{0,\tG} - (\tS w +R_h w, \partial_\tns v)_{0,\tG} \\
& \hskip 4cm + \gamma (h^{-1}(\tS w + R_h w), \tS v)_{0,\tG} \\
&= a_h(w,v) - (R_h w, \partial_\tns v)_{0,\tG} + \gamma (h^{-1} R_h w, \tS v)_{0,\tG} \\
& = a_h(v,w) + (\partial_\tns w, \nabla v \cdot \bs{d})_{0,\tG} - (\partial_\tns v, \nabla w \cdot \bs{d})_{0,\tG} \\
& \hskip 4cm - (R_h w, \partial_\tns v)_{0,\tG} + \gamma (h^{-1} R_h w, \tS v)_{0,\tG} \,,
\end{split}
\end{equation*}
where in the last step we have used Lemma \ref{lemma:non-symmetry}. We now pick $v = \psi := u-u_h$ and we use Lemma \ref{lemma:non-Gal-orth} with $v_h = w_I :=I_h w$ to get
\begin{equation}\label{eq:dua1}
\begin{split}
(\tilde{z},\psi)_{0,\tO} &= a_h(\psi, w) + E_1(\psi, w) + E_2(\psi,w)  \\
&= a_h(\psi, w- w_I) + E_1(\psi, w) + E_2(\psi,w) + E_3(u,w_I) \,,
\end{split}
\end{equation}
with
\begin{equation}\label{eq:def-E1}
E_1(\psi, w) := (\partial_\tns w, \nabla \psi \cdot \bs{d})_{0,\tG} - (\partial_\tns \psi, \nabla w \cdot \bs{d})_{0,\tG} \,,
\end{equation}
\begin{equation}\label{eq:def-E2}
E_2(\psi,w) := - (R_h w, \partial_\tns \psi)_{0,\tG} + \gamma (h^{-1} R_h w, \tS \psi)_{0,\tG}\,,
\end{equation}
\begin{equation}\label{eq:def-E3}
E_3(u,w_I) := (R_h u, \partial_\tns w_I)_{0,\tG} - \gamma (h^{-1}R_h u, \tS w_I)_{0,\tG}\,.
\end{equation}

We proceed to bound the four error terms on the right-hand side of \eqref{eq:dua1}. Concerning the first one, we observe that an argument similar to the one used in the proof of Proposition \ref{prop:estim-approx-err} yields the bound
\begin{equation}\label{eq:err-interp-w}
\Vert w-w_I \Vert_{W(\tO;\tT)} \leq C \, h_\Omega^{s_0} |\nabla w|_{\sigma_0,2,\Omega}\;.
\end{equation}
On the other hand, proceeding as in the proof of Proposition \ref{prop:continuity-a}, one gets
$$
|\, a_h(\psi,w-w_I)\, | \leq B \, \Vert \psi \Vert_{V(\tO;\tT)} \, \Vert w-w_I \Vert_{W(\tO;\tT)}
$$
for some constant $B>0$ independent of the meshsize. Hence, by  \eqref{eq:estim-discr-err-V} we deduce the existence of a constant $C_{0,0}>0$ independent of $u$, $\tilde{z}$ and the meshsize, such that
\begin{equation}\label{eq:first-bound-enh}
|\, a_h(\psi,w-w_I)\, | \leq C_{0,0} h_\Omega^{r+s_0} \, ( |\nabla u |_{s,p,\Omega} + K(u,f,g) ) \, \Vert \tilde{z} \Vert_{0,\tO} \,.
\end{equation}
Concerning the error term $E_1$ in \eqref{eq:dua1}, recalling \eqref{eq:d-smallness} we have
$$
|\,E_1(\psi, w) \, | \leq c_d \, h_\Gamma^{1/2+\zeta} \Vert \nabla w \Vert_{0,\tG} \Vert h^{1/2}\nabla \psi \Vert_{0,\tG}
$$
with $\Vert \nabla w \Vert_{0,\tG} \leq c (\Vert w \Vert_{1,\Omega} + |\nabla w|_{\rho,\Omega})$ by the trace theorem, and $\Vert h^{1/2}\nabla \psi \Vert_{0,\tG} \leq \bar{C}_I \, \Vert \psi \Vert_{V(\tO;\tT)}$ by \eqref{eq:scaled-5}. Hence, we deduce the existence of a constant $C_{0,1}>0$ independent of $u$, $\tilde{z}$ and the meshsize, such that
\begin{equation}\label{eq:second-bound-enh}
|\, E_1(\psi, w) \, | \leq C_{0,1} h_\Omega^{r + 1/2+\zeta}  \, ( |\nabla u |_{s,p,\Omega} + K(u,f,g) ) \, \Vert \tilde{z} \Vert_{0,\tO} \,.
\end{equation}
Next, considering the error term $E_2$ in \eqref{eq:dua1}, we have
\begin{equation*}
\begin{split}
|\,E_2(\psi, w) \, | &\leq  \Vert h^{-1/2} R_h w \Vert_{0,\tG} ( \Vert h^{1/2} \nabla \psi \Vert_{0,\tG} + \gamma \Vert h^{-1/2} \tS \psi \Vert_{0,\tG}) \\
& \leq c\, \Vert h^{-1/2} R_h w \Vert_{0,\tG} \Vert \psi \Vert_{V(\tO;\tT)}\,.
\end{split}
\end{equation*}
To estimate the norm involving $R_h w$, we can apply the equivalent statement of Lemma \ref{lem:consist-2} for $w$, with the summation index $p=2$. Thus,
denoting by $\sigma_0$ the constant defined in \eqref{eq:def-exponent-sigma} with $p=2$, we get
\begin{equation}\label{eq:bound-Rhw}
\Vert h^{-1/2} R_h w \Vert_{0,\tG}  \leq c \, h_\Gamma^{\sigma_0} \, ( |\nabla w |_{s_0,\Omega} + \Vert \tilde{z} \Vert_{0,\tO} )
\leq c' \, h_\Gamma^{\sigma_0} \, \Vert \tilde{z} \Vert_{0,\tO}
\end{equation}
Hence, we deduce the existence of a constant $C_{0,2}>0$ independent of $u$, $\tilde{z}$ and the meshsize, such that
\begin{equation}\label{eq:third-bound-enh}
|\, E_2(\psi, w) \, | \leq C_{0,2} h_\Omega^{r + \sigma_0}  \, ( |\nabla u |_{s,p,\Omega} + K(u,f,g) ) \, \Vert \tilde{z} \Vert_{0,\tO} \,.
\end{equation}

Finally, the error term $E_3$ in \eqref{eq:dua1} can be equivalently written as
\begin{equation*}
\begin{split}
E_3(u,w_I) & = (R_h u, \partial_\tns w)_{0,\tG} - \gamma (h^{-1}R_h u, \tS w)_{0,\tG} - (R_h u, \partial_\tns (w-w_I))_{0,\tG} \\
& \ \ \ +\gamma (h^{-1}R_h u, \tS (w-w_I))_{0,\tG} \,.
\end{split}
\end{equation*}
Recalling \eqref{eq:w-q}, one has
\begin{equation*}
\begin{split}
|\,E_3(u, w_I) \, |  &\leq  \Vert h^{-1/2} R_h u \Vert_{0,\tG} (h_\Gamma^{1/2} \Vert \nabla w \Vert_{0,\tG} +  \gamma \Vert h^{-1/2} R_h w \Vert_{0,\tG} + \\
& \qquad \qquad  \qquad \qquad + \Vert h^{1/2} \nabla (w-w_I) \Vert_{0,\tG} + \gamma \Vert h^{-1/2} \tS (w-w_I) \Vert_{0,\tG}) \\
& \hskip -1cm \leq \Vert h^{-1/2} R_h u \Vert_{0,\tG} (h_\Gamma^{1/2} \Vert \nabla w \Vert_{0,\tG} \!\! +  \gamma \Vert h^{-1/2} R_h w \Vert_{0,\tG} \!\! + c\, \Vert w - w_I \Vert_{W(\tO;\tT)}) \,.
\end{split}
\end{equation*}
Using \eqref{eq:est-consist-2}, \eqref{eq:bound-Rhw} and \eqref{eq:err-interp-w}, and observing that $\min(1/2, \sigma_0,s_0)=1/2$, we deduce the existence of a constant $C_{0,3}>0$ independent of $u$, $\tilde{z}$ and the meshsize, such that
\begin{equation}\label{eq:fourth-bound-enh}
|\, E_3(\psi, w) \, | \leq C_{0,3} h_\Omega^{\sigma+1/2 }  \, K(u,f,g)  \, \Vert \tilde{z} \Vert_{0,\tO} \,.
\end{equation}

We are ready to state our estimate on the $L^2$-norm of the discretization error $u-u_h$. Using the Riesz representation \eqref{eq:Riesz} together with \eqref{eq:dua1} and the bounds \eqref{eq:first-bound-enh}, \eqref{eq:second-bound-enh}, \eqref{eq:third-bound-enh}, and \eqref{eq:fourth-bound-enh}, we arrive at the following result.

\begin{thm}\label{theo:estim-L2-err}
Under the same assumptions for the validity of Theorem \ref{theo:estim-discr-err-a}, setting
$$
r_0:= \min(\, s_0, \, 1/2+\zeta, \, \sigma_0, \, \sigma-r+1/2 \,) \;,
$$
there exists a constant $C>0$ independent of $u$ and the meshsize such that
\begin{equation}\label{eq:estim-discr-err}
\Vert u - u_h \Vert_{0,\tO} \leq C \, h_\Omega^{r+r_0}\, ( |\nabla u |_{s,p,\Omega} + K(u,f,g) )\,.
\end{equation}
\end{thm}

\begin{rem}\label{rem:rate-L2}
{\rm
It is worth analyzing the predicted decay rate in different situations, making the exponent $r+r_0$ explicit.

In the most favorable case (smooth data and smooth boundary, hence, $\lambda_S=\infty$),
one has $r=s_0=1$, and $\sigma=3/2+2\zeta$, $\sigma_0=1+3\zeta/2$, whence $r+r_0= 3/2+\zeta$.

On the other hand, in the presence of boundary singularities (i.e., $1/2 < \lambda_S <1$), then $r$ and $s_0$ are any real numbers $< \lambda_S$, whereas $\sigma=\sigma_0=\lambda_S(1+\zeta)$, whence $r+r_0$ is any real number $< \lambda_S+\min(\lambda_S,1/2+\lambda_S\zeta)$.

Note that in both cases, the exponent $r+r_0$ is optimal ($=2$ or $<2\lambda_S$) provided $\zeta \geq 1/2$.
}
\end{rem}

\begin{rem}\label{re:non-optimality}
{\rm
Allowing $\zeta >0$ to be arbitrarily small is of paramount importance in complex engineering applications, as it essentially avoids restrictions in the mesh generation near the boundary.
In this case, the exponent in \eqref{eq:estim-discr-err} is sub-optimal.  Inspecting the behavior of the four error terms appearing in \eqref{eq:dua1} reveals that in the most favorable situation (smooth data and smooth boundary), the terms $a_h(\psi, w- w_I)$, $E_2(\psi,w)$, and  $E_3(u,w_I)$ exhibit the optimal decay $h_\Omega^2$ (see \eqref{eq:first-bound-enh}, \eqref{eq:third-bound-enh}, and \eqref{eq:fourth-bound-enh} with $p=2$). Thus, the responsible for sub-optimality is the term $E_1(\psi,w)$, which reflects the discrepancy between the normal $\bs{n}$ to the exact boundary $\Gamma$ and the normal $\tn$ to the surrogate boundary $\tG$.

To shed further light, assume that $\G$ is smooth and the mapping \eqref{eq:mappingM} is the closest-point projection, which implies $\bs{d}=\Vert \bs{d} \Vert \bs{n}$. Using the splitting $\partial_\tns = (\tn \cdot \bs{n}) \partial_n + (\tn \cdot \bs{t}) \partial_t$, after cancellation of the common term involving the product $\partial_n w \,\partial_n \psi$, one gets
$$
E_1(\psi,w) = \int_{\tG} \Vert \bs{d} \Vert \, (\tn \cdot \bs{t}) (\partial_n \psi \, \partial_t w - \partial_n w \, \partial_t \psi) \,.
$$
For sufficiently refined meshes, when $\tG$ tends to become parallel to $\G$, the factor $\tn \cdot \bs{t}$ may become small, and cancellations due to sign changes from one edge to the other along $\tG$ may occur. Indeed, the optimal rate $h_\Omega^2$ for $\Vert u - u_h \Vert_{0,\tO}$ is often observed in practical calculations with SBM.

Analogous considerations apply in the presence of geometric singularities.
}
\end{rem}

config: arxiv


TITLE: Evaluating $\int_{0}^{\pi} \frac{1}{\sqrt{u^2+2u\cos x +1}} \mathrm{d}x$
QUESTION [5 upvotes]: Does the integral $$\int_0^{\pi} \frac{1}{\sqrt{u^2+2u\cos x +1}} \text d x \hspace{30pt} (u \le 1) $$ has a closed form? If it has, how do we evaluate it?

I was solving a physics problem which I have asked on physics SE as well (here), and this integral popped out. How do I solve it?

REPLY [1 votes]: I prefer to write another answer for the case where $u>1$.
In the most general case,
$$\int_{0}^{\pi} \frac{1}{\sqrt{u^2+2u\cos x +1}} \,dx=\frac{2 }{u-1}K\left(-\frac{4 u}{(u-1)^2}\right)$$ and we can expand the rhs as a series built around $u=1$. This would give
$$\frac{2 }{u-1}K\left(-\frac{4 u}{(u-1)^2}\right)=\sum_{n=0}^\infty a_n \, u^n$$ where the first $a_n$'s are
$$\left\{d,\frac{1}{2}-\frac{d}{2},\frac{5 d}{16}-\frac{7}{16},\frac{17}{48}-\frac{7
   d}{32},\frac{169 d}{1024}-\frac{1775}{6144},\frac{14789}{61440}-\frac{269
   d}{2048},\frac{1781
   d}{16384}-\frac{100741}{491520},\frac{244835}{1376256}-\frac{3035
   d}{32768},\frac{338377
   d}{4194304}-\frac{276441853}{1761607680},\frac{1482803963}{10569646080}-\frac{59
   9569 d}{8388608}\right\}$$ where $d=\log \left(\frac{8}{u-1}\right)$.
Uisng the above terms, the result is quite decent up to $u=2$. For this value, the exact result is
$2 K(-8) \sim 1.68575$ while the expansion gives $1.68208$


TITLE: Integrating $\int \sec^2(x) \tan(x) dx$ by trig substitution
QUESTION [3 upvotes]: I know I am supposed to integrate
$$\int \sec^2(x) \tan(x) dx$$
by substituting $u = \tan(x)$ and get $du = \sec^2(x)$. However, why can't I use $u = \sec(x)$, $du = \tan(x) \sec(x)$?

REPLY [7 votes]: You can use either substitution.
If you use $u = \tan x$, $du = \sec^2 x \,dx$, then you get
$\displaystyle\int \sec^2 x \tan x \,dx = \int u\,du = \dfrac{1}{2}u^2+C = \dfrac{1}{2}\tan^2 x + C$ for some constant $C$.
If you use $u = \sec x$, $du = \sec x\tan x \,dx$, then you get
$\displaystyle\int \sec^2 x \tan x \,dx = \int u\,du = \dfrac{1}{2}u^2+C' = \dfrac{1}{2}\sec^2 x + C'$ for some constant $C'$.
These two answers are equivalent. Since $\sec^2 x = \tan^2 x + 1$, this just means that $C = C'+ \dfrac{1}{2}$.


TITLE: Satisfying Conditions for Inner Product in Space $R^2$
QUESTION [1 upvotes]: The solution states that this is an inner product iff $b=c$ and $d>b^2$. I am not sure how to solve/approach this problem. I know that there are 4 conditions for the inner product condition to satisfy. Can you explain how to use those four properties to find conditions for b,c, and d?
P.S. Properties of inner product

$<x,y>=<y,x>$
$<x+z,y>=<x,y>+<z,y>$
$<cx,y>=c<x,y>=<x,cy>$
$<x,x>=0$ if $x\neq0 $

REPLY [3 votes]: Eventually, you learn that positive definiteness will require that the determinant of the matrix of coefficients be positive. But you can see it directly by completing the square: With $x=(x_1,x_2)$,
$$\langle x,x \rangle = x_1^2 + 2bx_1x_2 + dx_2^2 = (x_1+bx_2)^2 + (d-b^2)x_2^2.$$
Can you finish now?


TITLE: differential of an integral (composition of functions)
QUESTION [0 upvotes]: I'm sorry, despite of all the hints here I don't manage to put it together.
Given $\frac{\mathrm{d}}{\mathrm{d}x}\int_{0}^{x}e^{(x-t)\,A}\,b(t)\,\mathrm{d}t$
As solution $b(x) + \int_{0}^{x}A\,e^{(x-t)\,A}\,b(t)\,\mathrm{d}t$ is sugged.
Apparently the chain rule applies somehow. But I don't know which functions are chained. Is it $\frac{\mathrm{d}}{\mathrm{d}x}\left[F(x)-F(0)\right]$ ? Doesn't seem right to me...

REPLY [0 votes]: Ah, I was terrified by the formula, but basically I just got to plug in my values:
Solution:
$\begin{align}f(x,x)\,\frac{\mathrm{d}}{\mathrm{d}x}(x) - f(x,0)\,\frac{\mathrm{d}}{\mathrm{d}x}(0) + \int_{0}^{x}\frac{\mathrm{d}}{\mathrm{d}x}f(x,t)\\[12pt]
 = e^{0}\,b(x)\,1-e^{x\,A}\,b(0)\,0 + \int_{0}^{x}A\,e^{(x-t)\,A}\,b(t)\,\mathrm{d}t \end{align}$


text:
\section{Introduction}
Complex engineering models are often described by systems of equations, where the model outputs depend on a set
of input parameters. Given values for the input parameters, the model output may be computed by solving a
discretized differential equation, which typically involves the solution (or series of solutions) of a system of linear
equations for the unknowns. Each of these computations may be very expensive depending on factors such as grid
resolution or physics components in the model. As the role of simulation gains prominence in areas such as
decision support, design optimization, and predictive science, understanding the effects of variability in the input
parameters on variability in the model output becomes more important. Exhaustive parameter studies (e.g.\ uncertainty
quantification, sensitivity analysis, model calibration) can be prohibitively expensive -- particularly for a large
number of input parameters -- and therefore accurate interpolation and robust surrogate models are essential.

We consider the model problem of a matrix equation whose matrix and right hand side depend on a set of
parameters. Let $s\in[-1,1]^d$ be a set of input parameters from an idealized parameter space. We equip the parameter
space with a bounded, separable, positive weight function $\omega:[-1,1]^d \mapsto \mathbb{R}_+$, where
$\omega(s)=\omega_1(s_1)\cdots\omega_d(s_d)$. (In a probabilistic context, this weight function represents a probability
measure on the input parameter space.) Let $A(s)$ be an $N\times N$ matrix-valued function where we assume that $A(s)$
is invertible for all $s\in[-1,1]^d$, and let $b(s)$ be the vector-valued function with $N$ square-integrable
components. We seek the vector valued function $x(s)$ that satisfies
\begin{equation}
\label{eq:parmmat}
A(s)x(s)=b(s),\qquad s\in[-1,1]^d.
\end{equation}
Such parameterized matrix problems often arise as an intermediate step when computing an approximate
solution of a complex model with multiple input parameters. They appear in such diverse fields
as differential equations with random (or parameterized) inputs~\cite{Babuska01,Prudhomme02}, electronic circuit
design~\cite{Li2009}, image deblurring models~\cite{Chung08}, and ranking methods for nodes in a
graph~\cite{Page99,Constantine10b}. Given a parameterized matrix equation, one may wish to approximate the vector
valued function that satisfies the equation or estimate its statistics. Once an approximation is constructed that is
cheaper to evaluate than the true solution, statistics -- such as mean and variance -- of the approximation represent
estimates of statistics of the true solution.

The approximation model is a vector of multivariate polynomials represented as a series of orthonormal
polynomial basis functions where each basis function is a product of univariate orthonormal polynomials. We employ the
standard multi-index notation; let $\alpha=(\alpha_1,\dots,\alpha_d)\in\mathbb{N}^d$ be a multi-index, and define
the basis polynomial
\begin{equation}
\pi_\alpha(s) = \pi_{\alpha_1}(s_1)\cdots\pi_{\alpha_d}(s_d).
\end{equation}
The polynomial $\pi_{\alpha_i}(s_i)$ is the orthonormal polynomial of degree $\alpha_i$, where the orthogonality is
defined with respect to the weight function $\omega_i(s_i)$. Then for $\alpha,\beta\in\mathbb{N}^d$,
\begin{equation}
\int_{[-1,1]^d} \pi_\alpha(s)\pi_\beta(s)\omega(s)\,ds =
\left\{
\begin{array}{cl}
1,&\alpha=\beta\\
0,&\mbox{otherwise}
\end{array}
\right.
\end{equation}
where equality between multi-indicies means component-wise equality. For a given index set $\sI\subset\mathbb{N}^d$ with
size $|\sI|<\infty$, the polynomial approximation can be written
\begin{equation}
\label{eq:polymodel}
x(s)\;\approx\;\sum_{\alpha\in\sI} \vx_\alpha \pi_\alpha(s) \;=\; \mX\bpi(s).
\end{equation}
The $N$-vector $\vx_\alpha$ is the coefficient of the series corresponding to $\pi_\alpha(s)$. The $N\times|\sI|$ matrix
$\mX$ has columns $\vx_\alpha$, and the parameterized vector $\bpi(s)$ contains the basis polynomials. The goal of the
approximation method is to compute the unknown coefficients $\mX$.

Such polynomial models have become popular for approximating the solution of PDEs with random inputs; they appear under
names such as polynomial chaos methods~\cite{Xiu02}, stochastic finite element methods~\cite{Ghanem91}, stochastic
Galerkin methods~\cite{Babuska04}, and stochastic collocation methods~\cite{Xiu05,Babuska05}. The Galerkin methods
compute the series coefficients such that the equation residual is orthogonal to the approximation space defined by the
index set $\sI$; they typically employ a full polynomial basis of order $n$ given by
\begin{equation}
\sI\;=\; \sI_n \;=\; \{\alpha\in\mathbb{N}^d\;:\;\alpha_1+\cdots+\alpha_d\leq n\},
\end{equation}
although such a basis set is not strictly necessary. The number of terms in this basis set is
$|\sI_n|=\binom{n+d}{n}$, which grows exponentially with the number of input parameters $d$. The process of computing
the coefficients involves solving a linear system of size $N|\sI|\times N|\sI|$, which can be prohibitively large for
even a moderate number of input parameters (6 to 10) and low order polynomials (degree $> 5$). For this reason, there
has been a flurry of recent work on solver strategies for the matrix equations arising from the Galerkin
methods~\cite{Pellissetti00,Rosseel08,Keese05,Elman07,Elman05,Ernst09}, including papers on
preconditioning~\cite{Powell09-2,Elman10,Powell08,Ullmann2010}. Such work has relied on knowing the matrix valued
coefficients $\mA_\alpha$ of a series expansion of the parameterized matrix $A(s)$,
\begin{equation}
\label{eq:Aexpansion}
A(s) = \sum_\alpha \mA_\alpha \pi_\alpha(s),
\end{equation}
which typically come from a specific form of the coefficients of the elliptic operator in the PDE models. Another
drawback of the Galerkin method is its limited ability to take advantage of existing solvers for the problem
$A(\lambda)x(\lambda)=b(\lambda)$ given a parameter point $\lambda\in[-1,1]^d$. In contrast, the pseudospectral and
collocation methods can compute coefficients of the model \eqref{eq:polymodel} using only evaluations of the solution
vector $x(\lambda)$. The chosen
points typically correspond to a multivariate quadrature rule, and computing the coefficients of the polynomial model
\eqref{eq:polymodel} becomes equivalent to approximating its Fourier coefficients with the quadrature
rule~\cite{Babuska05,Constantine09-2}. Thus, from the point of view of code reuse and rapid implementation, the
pseudospectral/collocation methods have a distinct advantage.

We propose a variant of the Galerkin method that alleviates the drawbacks of limited code reuse and memory
limitations. By formally replacing the integration in the Galerkin method by a multivariate quadrature rule -- a
step that is more often than not performed in practice -- we derive a decomposition of the linear system of equations
used to compute the Galerkin coefficients of \eqref{eq:polymodel}; such a method has been called Galerkin with
numerical integration in the context of numerical methods for PDEs~\cite{Canuto06}. This decomposition allows us to
compute the Galerkin coefficients using only evaluations of the parameterized matrix $A(\lambda)$ and parameterized
vector $b(\lambda)$ for points $\lambda\in[-1,1]^d$ corresponding to a quadrature rule. In fact, if one requires only
matrix-vector multiplies as in Krylov-based iterative methods (e.g.~\textsc{cg}~\cite{Hestenes52} or
\textsc{minres}~\cite{Paige75} solvers) for the Galerkin system of equations, this restriction can be relaxed to
matrix-vector products $A(\lambda)\vv$ for a given $N$-vector $\vv$; there is no need for the coefficients of an
expansion of $A(s)$ as in \eqref{eq:Aexpansion}. Therefore the method takes full advantage of sparsity of the
parameterized system resulting in reduced memory requirements. The decomposition also yields insights for
preconditioning the Galerkin system that generalize existing work. Additionally, the decomposition immediately reveals
bounds on the eigenvalues of the Galerkin system for the case of symmetric $A(s)$.

The remainder of the paper is structured as follows. In Section \ref{sec:galerkin}, we derive the Galerkin method for a
given problem and basis set. We then derive the decomposition of the matrix used to compute the Galerkin coefficients
and examine its consequences including bounds on the eigenvalues of the matrix, strategies for simple implementation
and code reuse, and insights into preconditioning. We then provide a numerical study of various preconditioners
suggested by the decomposition using the common test case of an elliptic PDE with parameterized coefficients in Section
\ref{sec:precon}; the codes for the numerical study are available in a \textsc{Matlab} suite of tools accessible online.
To emphasize the advantages of code reuse, we apply the method to an engineering test problem using existing codes in
Section \ref{sec:joe}. Finally, we conclude in Section \ref{sec:summary} with summarizing remarks.

\subsection{Notation}
For the remainder of the paper, we will use the bracket notation $\ip{\cdot}$ to denote a \emph{discrete} approximation
to the integral with respect to the weight function $\omega$. In other words, for functions
$f:[-1,1]^d\longrightarrow\mathbb{R}$,
\begin{equation}
\int_{[-1,1]^d}f(s)\omega(s)\,ds
\;\approx\;
\sum_{\beta\in\sJ} f(\lambda_\beta)\nu_\beta
\;\equiv\;
\ip{f},
\end{equation}
where the points $\lambda_\beta=(\lambda_{\beta_1},\dots,\lambda_{\beta_d})\in[-1,1]^d$ and weights
$\nu_\beta\in\mathbb{R}$ define a multivariate quadrature rule for multi-indicies in the set $\sJ\subset\mathbb{N}^d$.
We will abuse this notation by putting matrix and vector valued functions inside the brackets, as well. For example,
$\ip{A}$ denotes the mean of $A(s)$ computed with a quadrature rule.


TITLE: I need some help on this problem on Geometric Loci.
QUESTION [1 upvotes]: The base $BC$ of a variable triangle is fixed, and the sum $AB + BC$ is constant, the line drawn through the midpoint $D$ of $BC$ parallel to $AB$ meets the parallel $CP$ through $C$ to the internal bisector of the angle at $A$, in $P$. Prove that the locus of $P$ is a circle having $D$ for centre.
Here are my assumptions, since the base $BC$ if fixed, it follows that $AB$ must lie on a circle with centre $B$. But I don't know if this is gonna help as that would leave me with three possibilities for which $AB = BC$ or $AB < BC$ or $AB > BC$ and considering any of these cases seems to lead me to a dead end.

REPLY [2 votes]: My answer will assume $AB+AC$ is constant as opposed to the original statement (since by Jack's answer, under the original statement the locus of $P$ is not on a circle). Under this condition we can prove that the locus of $P$ is indeed a circle.
Let us find the length of $PD$ and try to establish it as a constant function of $BC$ and $AB+AC$. Let $E$ be the intersection of the angle bisector and line $PD$, and let $F$ be the intersection of the angle bisector and line $BC$.  Notice that $$\frac{PD}{CD} = \frac{DE}{DF} = \frac{AB}{BF}\implies PD = \frac{BC}{2}\cdot \frac{AB}{BF}$$
Because we also have$$\frac{AB}{AC}=\frac{BF}{FC}\implies \frac{AB}{AB+AC}=\frac{BF}{BF+FC}\implies \frac{AB}{BF} = \frac{AB+AC}{BC}$$
We know that $$PD = \frac{AB+AC}{2}$$which is constant. This implies that the locus of $P$ is a circle with center $D$ and radius $\frac{AB+AC}{2}$.


TITLE: Pressure=Heat? Else, Heat=Radiation?
QUESTION [3 upvotes]: A few days ago, I came up with this paradoxical conundrum:

Let's say we have a gas enclosed in a container. Supplying heat to it will cause an increase in the (magnitude of?) vibration of gas particles in the medium. Increasing the pressure of the gas in the container (say by decreasing its volume) will increase the extent of vibration too. So what differs pressure from heat?

After a bit of thought, I could find some potential flaws/counter-logics to my proposal -

Hotness is a sensation that a living body feels (by virtue of some receptor or something...)
Hence the physical expression of the supply of heat may be the same as that of the supply (excess) pressure.
The effect on the motion of gas molecules does not necessarily determine anything.

I think it'd be fair to admit though, the general "postulate" (if I am not wrong!?) that the supply of heat increases the randomness of the motion of particles in any medium.
This leads me to infer some interesting (but perhaps flawed) consequences -

When we feel hot, that sense must have something to do with the increase in the randomness of the motion of the "heat-receptor molecules" in the skin, and maybe their collision amongst themselves... (is this how they work?)

The way "heat receptor molecules" can be excited (vibrated), is by supply of radiant energy, for that would be the best way to get them excited fast (because otherwise your only option left is to mechanically vibrate them)

This means a patch of radiant energy is integral and necessary to feeling hot

One possible flaw in this silly little thought-experiment of mine could be that our heat receptors work chemically. Meaning that electron transfer plays the pivotal role in the sensation of heat, not inter-molecular collision.

Then again, some form of radiant energy is necessarily absorbed to push the electrons out of atoms and cause the biochemical reaction that occurs everytime we feel heat. This radiant energy can't be anything other than radiant heat.

The reason I propose is simple - if radiant heat is not the "electron-pusher", something else is, then heat is not a factor to feeling hot (as that something else may push the electron at any other time, even when heat isn't supplied, and we may feel hot!) Otherwise, let's consider heat is the electron pusher, but some other form of it - convection or conduction. Conduction and convection both rely on molecules moving. I guess molecules can't push electrons out! So it must be radiant heat, I believe...
As a final summary of my paradox, here are the riddles that are forced, if the aforesaid inferences be true -

If the electron-push theory of heat reception be true, touching a hot metal wouldn't feel anything - Heat is conducted through metals, with molecules vibrating and making the heat flow. Molecules can't vibrate and displace electrons, right?! So we must not feel heat at all!
Else if molecular collision theory of heat reception be true, pressure and heat would be the same to feel - Cause you know, both are, after all, molecules striking your skin! What's the difference?

Or else -

Every kind of heat - conduction or convection, or even radiation - principally flows as a "stream" of radiant energy flowing through a medium or without any medium at all. The side effects of this flow of radiant energy, when it is flowing through a medium, include molecular vibration (conduction), random molecular motion (convection) and electron-donation on supply of energy.

I used the word "stream" of radiant energy, because the "side effects" (as I proposed) seems to me almost as if small channels dividing out of a stream...
That basically summarises my proposal at the conclusion, that what essentially transfers heat through a metal rod, when held near a furnace, is never the vibration of the particles of the rod. But the flow of radiant heat of the furnace that goes through the rod.
My question is, obviously, am I right? How much? Where am I wrong?
P.S.- I feel I'm wrong in a lot of places, so just maybe enlist where I'm right; that'll reduce your workload...

REPLY [0 votes]: I am not an expert on the microscopic aspects, but I add the following.  Heat is associated with temperature and temperature is a measure of the average kinetic energy of the particles.  Pressure can do work and pressure is a measure of the average normal force per unit area of the particles on a surface.
A discussion of heat and work on the macroscopic level follows.
I think your example of the gas in the container needs clarification.  The result of heat or work (from pressure) is an increase in the internal energy of the gas. Heat and work are energy that is transferred across a system boundary; the gas has no heat and it has no work; it has internal energy that is changed by heat and or work done on the gas.
Here is a macroscopic discussion of the difference between Heat and Work. Heat and work are different, based on their definitions in thermodynamics. For a system, thermodynamics defines Heat as energy that crosses a system boundary- without mass transfer- solely due to a temperature difference. Work is energy that crosses a system boundary- without mass transfer- due to any intensive property difference other than temperature.  So work is a very broad concept, not just the force time distance concept of work used in mechanics.  [For example, mechanical work is a force acting through a distance; electrical work is electrical energy crossing the boundary of a system].   An "open" thermodynamic system is one with possible mass flow into/out of the system in addition to possible heat/work done on/by the system.
So pressure that moves a system boundary is associated with work and radiative heat transfer is associated with heat, and the effect of either on the system is a change in internal energy of the system.
For your example, if pressure from the surroundings acts quickly to compress the gas there is work done on the gas; there is little heat transfer between the surroundings and the gas for a rapid compression even if the surroundings and gas are at different initial temperatures (until later when the gas and the surroundings eventually reach the same temperature).  The temperature of the gas is increased by the change in internal energy of the gas from the work done on the gas. For an ideal gas, the internal energy is a function of temperature alone and by the ideal gas law an increase in pressure causes an increase in temperature.