TITLE: Inverse Laplace Transform with squared irreducible quadratic in denominator
QUESTION [0 upvotes]: Okay, the equation is $$\frac{2s^3-2s}{(4s^2-4s+5)^2}$$
So I use partial fractions with $$\frac{As+B}{4s^2-4s+5} + \frac{Cs+D}{(4s^2-4s+5)^2} = 2s^3-2s$$
and square that quadratic get $$(As+B)(16s^4-32s^3+56s^2-40s+25)+(Cs+D)(4s^2-4s+5) = 2s^3-2s$$
This is the part where I'm suppose to plug in values for s to solve for the variables. But the only one that would make part of it $=0$ would make all terms $=0$.
So I tried differentiating both sides some times, each time picking a easy value for s that would equate one variable to $0$. Then I end up with a system of equations to solve, but get the wrong answer. I don't know what I'm doing wrong. Would someone point out what I did wrong or try solving it for me and writing the system of equations they get?
Thanks!

REPLY [1 votes]: That should be $\frac{2s^3-2 s}{(4 s^2 - 4 s + 5)^2} = \frac{A s + B}{4 s^2 - 4 s + 5} + \frac{Cs + D}{(4 s^2 - 4 s + 5)^2}$.
Multiplying by the denominator, you should get
$ 2 s^3 - 2 s = (As + B) (4 s^2 - 4 s + 5) + C s + D $.
The "easy" values to substitute are the roots of $4 s^2 - 4 s + 5$, namely
$ 1/2 \pm i$.  Or you could use long division with remainder.


TITLE: Aren't all physical relations non-linear?
QUESTION [6 upvotes]: It is well-known that Hooke's Law is only approximately true and thus that linear relation
is merely an idealization not strictly corresponding to the reality. Wouldn't it be necessary/appropriate that all linear relations decribing physical phenomena be reformulated to contain non-linear terms for use in cases where higher accuracy is desired? In particular, what would oppose Maxwell's Equations eventually being modified to have non-linear terms? (Or put in other words, why these equations "must" be linear?)

REPLY [2 votes]: I think that all physics is linear (or some other order) approximation to laws of nature. What I mean is that all formulas that are linear (or some other order) like gravity F = M*m*G/(R^2) in fact should be non linear functions like F = F(M,m,R). The only reason why it appears not so is that we perceive world and do experiments in very limited range of orders of magnitude of variables. Probably if you increase R much beyond or below usual observation you will see that gravity is not proportional to R^(-2). If gravity works differently on the scale R similar to size of universe then you won't detect that at all.
"In particular, the action of any observable (as well as any evolution operator or any other physically meaningful transformation) on the wave function is always exactly linear."
That's true in the theory that was based on limited experiments. Experiments performed for different scales or energies will contradict the theory and open doors for some new theory.

REPLY [1 votes]: Let me note that (under some conditions) a nonlinear theory can be precisely (rather than approximately) described by a linear theory using, e.g., Carleman linearization (embedding). See, e.g., http://engineering.ucsb.edu/~bamieh/pubs/doylefest_seq.pdf , p.38). The price you have to pay - the linear theory "lives" in an infinite-dimensional space.


TITLE: Fourier expansion of Takagi-function (everywhere non differentiable function).
QUESTION [7 upvotes]: Let us consider Takagi-function defined by
$T(x) \colon\!= \sum_{n=0}^{\infty}s(2^nx)/2^n$,
where $s(x) \colon\!\!= \underset{n \in {\Bbb Z}}{\mathrm{min}} \,|x-n|$.
$T(x)$ has its period $1$, so I am trying to develop its Fourier expansion
$S(T)(x) \colon\!=\sum_{n= -\infty}^{\infty}\left(a_n\,sin(2 \pi n x) +
b_n\,cos(2 \pi n x)\right)$.
Q: Where do we have $T(x) = S(T)(x)$?
For the case of Weierstrass function $W(x) \colon\!= \sum_{n=0}^{\infty}a^n cos(b^n \pi x)$ with $ab > 1 + 3\pi/2$, $ 0 < a < 1$ and $b$ an odd integer, it is already expanded in Fourier series with period $2$ by shape.
By Carlson, we know $T(x) = S(T)(x)$ except for sets of Lebesgue measure $0$. I would like to know whether $T(x) = S(T)(x)$ holds everywhere for Takagi-function as well.

REPLY [5 votes]: Yes. The triangle wave $s(x):=\min_{k\in\mathbb{Z}}\big|x-k\big|$ has an absolutely convergent Fourier series
$$s(x)=\frac{1}{4}-\frac{2}{\pi^2}\sum_{k=0}^\infty\frac{1}{(2k+1)^2}\cos\big(2\pi (2k+1)x\big)\, ,\qquad x\in\mathbb{R}.$$
Therefore
$$T(x):=\sum_{n=0}^\infty2^{-n}s(2^nx)= \frac{1}{4}\sum_{n=0}^\infty2^{-n}-\frac{2}{\pi^2}\sum_{n=0}^\infty\sum_{k=0}^\infty \frac{1}{2^n(2k+1)^2}\cos\big(2\pi 2^n(2k+1)x\big)\, .$$
By absolute convergence we can reorder the double sum. Since every positive integer $m\ge1$ writes uniquely as $m=2^n(2k+1)$ for nonnegative integers $n$ and $k$, we then obtain an absolutely convergent  Fourier series for $T(x)$:
$$T(x)=\sum_{m=0}^\infty a_m\cos(2\pi m x)$$
with $a_0=1/2$ and for $m\ge 1$
$$a_m:=-\frac{2^{\nu(m)+1}}{\pi^2m^2}\, , $$
where $2^{\nu(m)}$ is the maximum power of $2$ that divides $m$.


TITLE: How does $\sum^n_{i=1}\sum^n_{j=1}x_ix_j\mathrm{Cov}\left(Z_i,Z_j\right)$ equal to $\mathrm{Cov}\left(\sum^n_{i=1}x_iZ_i, \sum^n_{j=1}x_jZ_j\right)$?
QUESTION [1 upvotes]: I have a proof which I don't understand quite one step:
Let $\sum$ be a covariance matrix. Then $\sum$ is non-negative definite.
Proof:
Let $Y=x_1Z_1+x_2Z_2+...+x_nZ_n$. We want to prove $\underline{x}^T\sum\underline{x} \geq0$ for all $\underline{x} \in \Bbb R^n\backslash\{0\}$
Then
$$\underline{x}^T\sum\underline{x}=\sum^n_{i=1}\sum^n_{j=1}x_ix_j\operatorname{Cov}(Z_i,Z_j)=\operatorname{Cov}(\sum^n_{i=1}x_iZ_i, \sum^n_{j=1}x_jZ_j)=\operatorname{Cov}(Y,Y)=\operatorname{ Var}(Y)\geq0$$
Therefore, $\sum$ is n.n.d.
I don't quite understand the part where $$\sum^n_{i=1}\sum^n_{j=1}x_ix_j\operatorname{Cov}\left(Z_i,Z_j\right)=\operatorname{Cov}\left(\sum^n_{i=1}x_iZ_i, \sum^n_{j=1}x_jZ_j\right)$$
I know the formula $\operatorname{Cov}(au+bv,w)=a\operatorname{Cov}(u,w)+b\operatorname{Cov}(v,w)$
So, to understand how how they derived it, I tried to derive it myself by expanding $\sum^n_{i=1}\sum^n_{j=1}\operatorname{Cov}\left(Z_i,Z_j\right)$:
$$x_1x_1\operatorname{Cov}(Z_1,Z_1)+x_1x_2\operatorname{Cov}(Z_1,Z_2)+x_1x_3\operatorname{Cov}(Z_1,Z_3)+...+x_1x_n\operatorname{Cov}(Z_1,Z_n)+$$
$$+x_2x_1\operatorname{Cov}(Z_2,Z_1)+x_2x_2\operatorname{Cov}(Z_2,Z_2)+x_2x_3\operatorname{Cov}(Z_2,Z_3)+...+x_2x_n\operatorname{Cov}(Z_2,Z_n)+$$
$$+$$
$$.$$
$$.$$
$$+$$
$$+x_nx_1\operatorname{Cov}(Z_n,Z_1)+...+x_nx_n\operatorname{Cov}(Z_n,Z_n)$$
So, $$\underbrace{x_1x_1}_{a}\operatorname{Cov}(\underbrace{Z_1}_{w},\underbrace{Z_1}_{u})+\underbrace{x_1x_2}_{b}\operatorname{Cov}(\underbrace{Z_1}_{w},\underbrace{Z_2}_{v})=\operatorname{Cov}(x_1x_2Z_1+x_1x_2Z_2,Z_1)$$
Ok so I iterate:
$$\underbrace{1}_{a}\operatorname{Cov}(\underbrace{x_1x_2Z_1+x_1x_2Z_2}_{u},\underbrace{Z_1}_{w})+\underbrace{x_1x_3}_{b}\operatorname{Cov}(\underbrace{Z_1}_{w},\underbrace{Z_3}_{v})=\operatorname{Cov}(x_1x_1Z_1+x_1x_2Z_2+x_1x_3Z_3,Z_1)$$
$$.$$
$$.$$
$$\operatorname{Cov}(x_1(\sum^n_{i=1}x_iZ_i),Z_1)$$
So the whole sum becomes $$\sum^n_{i=1}\operatorname{Cov}\left(x_i\left(\sum^n_{j=1}x_jZ_j\right),Z_i\right)$$
So how does this equal to $$\operatorname{Cov}\left(\sum^n_{i=1}x_iZ_i, \sum^n_{j=1}x_jZ_j\right)?$$

REPLY [1 votes]: Here is a short version of this proof:
\begin{align}
x^{\top}\Sigma x&=\mathsf{E}[x^{\top}(Z-\mathsf{E}Z)(Z-\mathsf{E}Z)^{\top}x]  \\
&=\mathsf{E}[Z^{\top}x-\mathsf{E}Z^{\top}x]^2=\operatorname{Var}(Z^{\top}x).
\end{align}


TITLE: The units digit of $1!+2!+3!+4!!+5!!+\dots+k\underset{\left \lfloor \sqrt{k} \right \rfloor \text{ times}}{\underbrace{!!!\dots!}}$
QUESTION [10 upvotes]: For natural numbers $n\ge m$, let $n\underset{m \text{ times}}{\underbrace{!!!\dots!}}=n(n-m)(n-2m)(n-3m)\dots$ where all factors are natural numbers (we exclude $0$ and negative factors).
Question:

What is the units digit of $1!+2!+3!+4!!+5!!+\dots+k\underset{\left
\lfloor \sqrt{k} \right \rfloor \text{
times}}{\underbrace{!!!\dots!}}+\dots+1992\underset{44 \text{
times}}{\underbrace{!!!\dots!}}$? ($\left \lfloor \cdot \right
\rfloor$ denotes the floor funtion).

My Attempt (Is wrong as  Peter Foreman commented below):
Consider the first $9$ terms:
$1!+2!+3!+4!!+5!!+6!!+7!!+8!!+9!!!$
$=1+2+6+8+15+48+105+384+162=731$
Each of the remaining terms includes at least on factor that ends with $0$. Therefore, the each term ends with $0$.
Hence the units digit of the given expression is equal to the units digit of the sum of the first $9$ terms. So, $1$ is the units digit of the given expression.

Peter Foreman said: "$17!!!!=9945$". This showed me that my attempt is wrong. Thanks  Peter Foreman.

Any help would be appreciated. THANKS.

REPLY [2 votes]: When $k \geq 25$ and the floor part $p$ of $k^{0.5}$ is coprime with $10$, the unit digit of $k’=k! \ldots !$ is $0$ (there is an even number in $k,k-p$ and one divisible by $5$ in $k!!!!!$).
When $k \geq 25$ and $p \wedge 10=2$, there is going to be a number divisible by $5$ in $k,k-p,k-2p,k-3p,k-4p$, and $k’$ is divisible by $5$ and congruent to $k$ mod $2$, so the unit digit of $k’$ is $5$ if $k$ is odd and $0$ if $k$ is even.
When $k \geq 25$ and $p \wedge 10=5$, then $k(k-p)$ is even, so $k’$ is even. The congruence mod $5$ is trickier: $k’$ is congruent to $k^r$ mod $5$, where $r$ is the number of factors in he product (ie $r-1$ is the floor part of $(k-1)/p$, so $r$ is either $p$, $p+1$ or $p+2$).
When $10 | p$, as above, $k’$ is congruent to $k^r$ mod $10$.
Note that everything depends only on the unit digit of $k$, $p$ and $r$: when $p$ is set, and $p \wedge 10=2$, the sum of any four $k’$ corresponding to consecutive $k \geq 25$ vanishes mod $10$.
When $p,r$ are set and $p \wedge 10=5$, the sum of any five $k’$ corresponding to consecutive $k \geq 25$ is always divisible by $10$.
When $p,r$ are set and $p \wedge 10=10$, the sum of any ten $k’$ corresponding to consecutive $k \geq 100$ is congruent to $3$ mod $10$ if $4|r$ and $5$ mod $10$ otherwise.
So all we need now is time for processing all integers from $1$ to $1992$.


TITLE: Using mathimagical induction to prove Inverse rule for groups elements
QUESTION [4 upvotes]: The question is:

Use mathematical induction to prove that if  $a_1,a_2,\ldots,a_n$ are elements of a group $G$ then $$(a_1a_2\cdots a_n)^{-1} = a_n^{-1}a_{n-1}^{-1}\cdots a_1^{-1}$$

May attempts:
Starting with a base case where $k=2$ we have $a_1$,$a_2~~$ $(a_1,a_2)^{-1}=e$
we multiply out by $(a_2)^{-1} (a_1)^{-1} a_1,a_2 (a_1,a_2)^{-1}= (a_2)^{-1} (a_1)^{-1} e$
and get
$e*e*(a_1,a_2)^{-1}=(a_2)^{-1}(a_1)^{-1}$
now i assume true for the $k$-th case and want to show it for the $k+1$ case
${(a_1,a_2,\ldots,a_k)^{-1}} = (a_k)^{-1},\ldots,(a_1)^{-1}$
now i multiply by $(a_{k+1})^{-1}$
$(a_{k+1})^{-1} * {(a_1,a_2,\ldots,a_k)^{-1}} = (a_{k+1})^{-1} * (a_k)^{-1},\ldots,(a_1)^{-1}$
My main question being is can i treat $(a_{k+1})^{-1}$ as one item and ${(a_1,a_2,\ldots,a_k)^{-1}}$ as anther item say $(a_d)^{-1}$ and convert it into  $(a_d,a_{k+1} )^{-1} = (a_{k+1})^{-1} * (a_k)^{-1},\ldots,(a_1)^{-1}$
Does this finish the proof?

REPLY [1 votes]: Hint $\ $ Below is the essence of the inductive step, with $\rm\: x' = x^{-1}$
$\rm\begin{eqnarray}\qquad  &&\ \rm (ab\cdots\! &\rm y& &\rm\!\! z)  (z'& &\rm\! y'&\rm\! \cdots b'a')\\
\rm &=&\ \rm  (ab\cdots\! &\rm y&\!\!)&\rm\! z\ \ z'&\!\!(&\rm\! y'&\rm\!\cdots b'a')\quad by\  associativity\\
\rm &=&\ \rm  (ab\cdots\! &\rm y&\!\!)&\ 1&\!\!(&\rm\! y'&\rm\!\cdots b'a')\quad by\, \ zz' = 1\\
\end{eqnarray}$
The prior $ = 1$ by the induction hypothesis, after absorbing the $1$. Ditto in the other order.


TITLE: Show that the $n\times n$ matrices with determinant $=1$ forms a $C^1$ surface of dimension $n^2-1\in \mathbb{R^{n^2}}$
QUESTION [3 upvotes]: I am told that I need to find a path $c(t)$ such that $c(t)=x(t), X(0)=x \forall X s.t.  det X=1$. So I can show that $d/dt(f(c(t))$ at $t=0=[d_{f(c(t))}f](c'(t))]\ne 0$
My problem is how to explicitly construct such a path $c(t)$?

REPLY [0 votes]: If I understand your question well, what you are looking for is the following: Given a matrix $X$ with $\det(X)=1$ find an explicit curve $c$ with $c(0)=X$ such that $t\mapsto \det(c(t))$ has non-vanishing derivative for $t=0$ (so this shows that $D\det(X)\neq 0$. If this is indeed the question, then you can simply take $c(t)=(1+t)X$, since this has determinant $(1+t)^n$, so its derivative at $t=0$ is $n\neq 0$.


TITLE: Lower bound on the minimum eigenvalue of sum of two matrices
QUESTION [2 upvotes]: Assume that $A$ is a symmetric positive definite matrix and $B$ is a symmetric (can potentially negative entries). Is the following bound correct?
$$\lambda_{min}(A+B)\geq \lambda_{min}(A)+\lambda_{min}(B)$$

REPLY [1 votes]: The bound also follows from an application of Weyl's inequalities (sometimes called the Courant-Weyl inequalities).


TITLE: Alpha particles and quantum tunneling
QUESTION [0 upvotes]: collapse and revival
There are a few things that I quite didn't understand about the tunneling of alpha particles.

Where does the kinetic energy of the alpha particle comes from? Is it because of the electromagnetic repulsion between protons?

There's a small probability of finding the particle inside the Coulomb barrier. What would that mean? Would the particle return to the nucleus or repelled away by the electromagnetic force?

REPLY [0 votes]: Where does the kinetic energy of the alpha particle comes from? Is it because of the electromagnetic repulsion between protons?

Kinetic and potential energies do not commute - hence, according to the uncertainty principle, the particle always have some kinetic energy. E.g., a particle bound in an infinite square well has only kinetic energy - different at every energy level.
However, even classically, being in a bound state does not mean zero kinetic energy - e.g., a planet rotating around the Sun has kinetic energy, yet not sufficiently high to escape. Also, as per classical physics, kinetic energy increases as the particle moves to the regions with lower potential energy.

There's a small probability of finding the particle inside the coulomb barrier. What would that mean? Would the particle return to the nucleus or repelled away by the electromagnetic force?

It means that, when measuring the position of the particle, one could also get values where it is inside the barrier. We consider a particle as a wave packet - i.e., a superposition of the Hamiltonian eigenstates - which is initially localized within the nucleus. As the system evolves, at times $t>0$ the superposition is no more localized and we have finite probability of finding particle elsewhere, including inside the barrier and outside the nucleus. In principle, since every eigenstate evolves periodically as $e^{-iEt/\hbar}$, one could expect that after some time all the states return to the situation that we had at $t=0$ - this is known as collapse and revival of the wave function. However, incommensurability of the periods makes the revival time very long... and when we take continuum limit for energy states, it becomes infinite - so particle never really returns to the initial state (although we always have non-zero probability of finding it in the nucleus).


TITLE: Derivative using chain rule with multi variable function
QUESTION [0 upvotes]: How would I calculate the total derivative using the chain rule with a multi-variable function? I need a total derivative with respect to $(t)$
$$R(t) = \max_x f(x^*(t, p), t)$$
Is this correct?
$$
\frac{\partial R}{\partial t} = \frac{\partial f}{\partial x} \frac{\partial x}{\partial t} \frac{\partial x}{\partial p} + \frac{\partial f}{\partial t}
$$
Or do I need to extend it out?
$$
\frac{\partial R}{\partial t} = \frac{\partial f}{\partial x} \frac{\partial x}{\partial t} + \frac{\partial f}{\partial x}\frac{\partial p}{\partial t} + \frac{\partial f}{\partial t}
$$

REPLY [1 votes]: I'm not sure how you got to your first option, but your second option is almost right. The only thing missing is the partial of $x^*$ with respect to $p$ in the second term.
\begin{align}
\frac{\mathrm d R}{\mathrm d t} &= \frac{\partial f}{\partial x^*} \frac{\mathrm d x^*}{\mathrm d t} + \frac{\partial f}{\partial t}\\
&= \frac{\partial f}{\partial x^*} \left( \frac{\partial x^*}{\partial t} + \frac{\partial x^*}{\partial p}\frac{\partial p}{\partial t} \right) + \frac{\partial f}{\partial t}\\
&= \frac{\partial f}{\partial x^*} \frac{\partial x^*}{\partial t} + \frac{\partial f}{\partial x^*}\frac{\partial x^*}{\partial p}\frac{\partial p}{\partial t} + \frac{\partial f}{\partial t}
\end{align}


TITLE: Error Estimation and Propagation through Trigonometric Functions
QUESTION [1 upvotes]: I have a problem whereby I have to estimate the absolute error of a value calculated using values read from sensors. The equation used to calculate the value is $T_{xy} = R \tan \theta$.
The sensor that provides the $R$ value has a fixed absolute error value $\epsilon_{R}=\pm 5$ m. However, the sensor that provides $\theta$ values only provides a standard deviation value, $\sigma_{\theta}$. The first thing that I need to do is use this to estimate $\epsilon_{\theta}$, the absolute error in $\theta$.
My logic for this estimate is to use $\epsilon_{\theta} = 3\sigma_{\theta}$, as if the value of $\theta$ is normally distributed $3\sigma_{\theta}$ will cover 99% of the probably values. I'm not sure if this is a good estimate to use, I'd welcome suggestions on this.
Now that I have my values $\epsilon_{R}$ and $\epsilon_{\theta}$ I need to estimate $\epsilon_{T_{xy}}$. To do this I think I need to use the formula $\epsilon_{T_{xy}}=T_{xy}(\frac{\epsilon_{R}}{R} + \frac{\epsilon_{\tan\theta}}{\tan\theta})$. Is this correct?
To calculate the above formula I need to define $\epsilon_{\tan\theta}$. I have defined this as $\epsilon_{\tan\theta} = |\sec^{2}\theta|\epsilon_{\theta}$. I am unsure whether this is correct.
It would be great if somebody could either validate my approach, or point out any errors, as I'm currently getting error estimates that make my device unusable. The dominant component of the estimate seems to be the error associated with the angle $\theta$.

REPLY [1 votes]: What you need is to figure out the distribution of $T_{xy}$ given distributions for $R$ and $\tan \theta$. I'll assume that $R$ and $\theta$ are Gaussian variables with standard deviations given by the stated errors or maybe half of the stated errors (you would have to try to puzzle it out from whatever product literature you have). Note that this is already demonstrably incorrect, since $\theta$ has a finite range and a Gaussian variable must have an infinite range.
There are at least three ways to go from there:
(1) Approximate $R \tan \theta$ as linear in $R$ and $\theta$. This is inaccurate because you're ignoring higher-order terms, and it probably becomes more inaccurate as $\theta$ is farther from 0. However, it does make the problem much simpler, and it may be the only way to get an exact formula. A linear combination of Gaussian variables is again Gaussian. (Even if you don't assume that $R$ and $\theta$ are Gaussian, assuming linearity is still a big simplification.) The mean of a linear combination of Gaussian variables is a linear combination of the means and the variance is a linear combination of the variances.
(2) Try harder to get an exact formula. Use the so-called change of variables method. That isn't guaranteed to work.
(3) Use a Monte Carlo method to approximate the distribution of $R \tan \theta$ with a histogram.
My advice is to try the Monte Carlo method first, in order to get a feeling for the problem, and then the linearization method. I can help you work out the formulas for the linearization if you want. If the linearization method is too inaccurate, consider using the Monte Carlo method -- does it matter how long it takes to compute the error estimate? If not, then generating a million numbers via Monte Carlo won't hold you back.
EDIT: Here's a derivation of a linearized approximation and Monte Carlo approximation for comparison. In summary, mean and standard deviation of the linearized Gaussian approximation are:
$$\mu_{T_{x,y}}=R_{0}\,\tan \theta_{0}$$
$$\sigma_{T_{x,y}}=\sqrt{\tan ^2\theta_{0}\,\sigma_{R}^2+R_{0}^2\,
 \tan ^4\theta_{0}\,\sigma_{\theta}^2+2\,R_{0}^2\,\tan ^2
 \theta_{0}\,\sigma_{\theta}^2+R_{0}^2\,\sigma_{\theta}^2}$$
Maxima session:
(%i2) T[x,y]:R*tan(theta)
(%i3) foo:taylor(T[x,y],[R,theta],[R_0,theta_0],[1,1])
(%o3)/T/ tan(theta_0) R_0 + (tan(theta_0) (R - R_0)
                                 2
                           + (tan (theta_0) + 1) R_0 (theta - theta_0)) + . . .
(%i4) foo:subst([R_0 = R[0],theta_0 = theta[0]],expand(foo))
                            2                                         2
(%o4) tan(theta ) R + R  tan (theta ) theta + R  theta - theta  R  tan (theta )
               0       0           0           0              0  0           0
                                                                    - theta  R
                                                                           0  0
(%i5) A:coeff(foo,R)
(%o5)                             tan(theta )
                                           0
(%i6) B:coeff(foo,theta)
                                   2
(%o6)                        R  tan (theta ) + R
                              0           0     0
(%i7) C:expand(-B*theta-A*R+foo)
                                    2
(%o7)                - theta  R  tan (theta ) - theta  R
                            0  0           0         0  0
(%i8) mu[T[x,y]]:expand(C+B*theta[0]+A*R[0])
(%o8)                           R  tan(theta )
                                 0          0
(%i9) sigma[T[x,y]]:sqrt(expand(B^2*sigma[theta]^2+A^2*sigma[R]^2))
              2              2    2    4              2
(%o9) sqrt(tan (theta ) sigma  + R  tan (theta ) sigma
                     0       R    0           0       theta
                                     2    2              2        2      2
                                + 2 R  tan (theta ) sigma      + R  sigma     )
                                     0           0       theta    0      theta
(%i10) [R[0],sigma[R],theta[0],sigma[theta]]:[17.29,0.25,0.9,0.01]
(%i11) ev([A,B,C])
(%o11)    [1.260158217550339, 44.7464980980593, - 40.27184828825338]
(%i12) load(distrib)
(%i13) load(descriptive)
(%i14) R_sample:random_normal(R[0],sigma[R],10000)
(%i15) theta_sample:random_normal(theta[0],sigma[theta],10000)
(%i16) T_xy_sample:R_sample*tan(theta_sample)
(%i17) histogram(T_xy_sample,nclasses = 20)
(%i18) [mean(T_xy_sample),ev(mu[T[x,y]])]
(%o18)              [21.79822439822625, 21.78813558144536]
(%i19) [std(T_xy_sample),ev(sigma[T[x,y]])]
(%o19)             [0.5516597186615362, 0.5472429351144614]
(%i20) tex('(mu[T[x,y]]) = mu[T[x,y]])
$$\mu_{T_{x,y}}=R_{0}\,\tan \vartheta_{0}$$
(%o20)                               false
(%i21) tex('(sigma[T[x,y]]) = sigma[T[x,y]])
$$\sigma_{T_{x,y}}=\sqrt{\tan ^2\vartheta_{0}\,\sigma_{R}^2+R_{0}^2\,
     \tan ^4\vartheta_{0}\,\sigma_{\vartheta}^2+2\,R_{0}^2\,\tan ^2
     \vartheta_{0}\,\sigma_{\vartheta}^2+R_{0}^2\,\sigma_{\vartheta}^2}$$

Maxima input script:
T[x, y] : R*tan(theta) $
    foo : taylor (T[x, y], [R, theta], [R_0, theta_0], [1, 1]);
    foo : subst ([R_0 = R[0], theta_0 = theta[0]], expand (foo));
    A : coeff (foo, R);
    B : coeff (foo, theta);
    C : expand (foo - A*R - B*theta);
    mu[T[x, y]] : expand (A*R[0] + B*theta[0] + C);
    sigma[T[x, y]] : sqrt (expand (A^2*sigma[R]^2 + B^2*sigma[theta]^2));
    [R[0], sigma[R], theta[0], sigma[theta]] : [17.29, 0.25, 0.9, 0.01] $
ev ([A, B, C]);
load(distrib) $
    load(descriptive) $
R_sample : random_normal (R[0], sigma[R], 10000) $
    theta_sample : random_normal (theta[0], sigma[theta], 10000) $
T_xy_sample : R_sample * tan(theta_sample) $
    histogram (T_xy_sample, nclasses = 20) $
[mean(T_xy_sample), ev (mu[T[x,y]])];
[std(T_xy_sample), ev (sigma[T[x,y]])];
tex ('(mu[T[x, y]]) = mu[T[x, y]]);
tex ('(sigma[T[x, y]]) = sigma[T[x, y]]);


TITLE: Solving Difference Equation
QUESTION [0 upvotes]: I have the following equation:
$$\frac{Z-1}{{Z^2}-2Z+6}$$
I tried to apply difference equation on it:
$$y_{n+2}-5y_{n-1}-6y_n = x_{n+1}-x_n$$
I need to find equation expressed as $y_n$
I have the answer (At the end of the chapter) to be:
$$y_n = 5y_{n-1}-6y_{n-2}+x_{n-1}-x_{n-2}$$
Can somone please tell me how did in this case the second equation expressd as $y_n$ ? I cant figure it out!

REPLY [1 votes]: The differences between your equation and the book's are:
the index has been shifted by 2
the sign on the 6 has been corrected
two terms have been moved to the other side
To see the index shift, define $m=n+2$ and substitute it in.
Also, the $2$ in the denominator has become $5$ in both your solution and the book's.  Did you copy it correctly?


TITLE: compact set boundary condition prove
QUESTION [1 upvotes]: the question is:

Prove that a set E ⊂ R is open if and only if E is not contains any of its boundary points.

I'm kind of not sure with good enough(?) proof... so please look at my several approach and fixed it. and if you have your own proof please let me know

E is an open but it does not contains all its boundary points.
E is open so define the E as
                 E = (a,b]
at point b it as the b+ϵ, so it does not contains all its boundary points.
so  'E is an open but it does not contains all its boundary points. ' assumption is      contradicted

REPLY [0 votes]: $BdE=\overline E-intE=\overline E \cap \overline {(\Bbb R-E)}$ and $BdE\cap intE=\emptyset$. If $E$ is open then $intE=E$ and thus $BdE\cap E=\emptyset$.
Now if $BdE\cap E=\emptyset$ then $E\cap \overline E \cap \overline {(\Bbb R-E)}= E \cap \overline {(\Bbb R-E)}=\emptyset$. Thus $E\subset \Bbb R-\overline {(\Bbb R-E)}=intE$ and so $E=int E$


TITLE: Prove that f is a constant function.
QUESTION [0 upvotes]: My friend asked me for help on the following question:
If $f:\mathbb R \to \mathbb Z$ is a continuous function, prove that $f$ is a constant function.
I tried the following proof, but I am generally bad at writing proofs, so I would like to know if the following is acceptable:
Proof
Suppose $f$ isn't constant, then $\exists x_1, x_2 \in \mathbb R$ such that $x_1 \ne x_2$ and $f(x_1) \ne f(x_2)$.
Note that $f(x_1)$, $f(x_2) \in \mathbb Z$.
Assume, without loss of generality, that $f(x_2) > f(x_1)$.
Then $f(x_2) - f(x_1) \ge 1$ and so $f(x_2) > f(x_2) - 0.5 > f(x_1)$.
But $f(x_2) - 0.5 \not\in \mathbb Z$, so $\not \exists x \in \mathbb R$ such that $f(x) = f(x_2) - 0.5$.
By the Intermediate Value Theorem this means that $f$ isn't continuous; a contradiction.
Therefore $f$ is constant.

REPLY [1 votes]: One version of the Intermediate Value Theorem just says the image of a connected subset of $\mathbf R$ is connected. Now in a discrete topological space, the only connected subsets are points. Hence the image of $\mathbf R$ is a point.


TITLE: Linear Algebra - Find the line that passes through a point and intersects two lines.
QUESTION [2 upvotes]: Basically I have to find a line that passes through a point and intersects two lines.
The Point $P=(6,-2,-2)$
The first Line is
$$x = 1+t$$
$$y = 5 + 2t$$
$$z = 3 +4t$$
The second line is
$$x = 3 + 3s $$
$$y = 1 + s $$
$$z = 4 + 2s $$
I have to express the answer in parametric form:
$$x = 6 + Aw $$
$$y = -2 + Bw $$
$$z = -2 + Cw $$
I basically have to find $A$, $B$ and $C$.
Thanks!

REPLY [1 votes]: The plane containing all lines through $P$ and some point on the first line has normal vector
$$
\vec n_{\alpha}=
\begin{pmatrix}
1-6\\
5-(-2)\\
3-(-2)
\end{pmatrix}
\times
\begin{pmatrix}
1\\
2\\
4
\end{pmatrix}
$$
When this cross product has been computed, a similar computation can be carried out to determine a normal vector $\vec n_{\beta}$ for the plane containing the point and the second line. Then a direction vector for the line of intersection can be found as
$$
\vec r=\vec n_{\alpha}\times\vec n_{\beta}
$$
and the coordinates of $\vec r$ should contain the figures $A,B,C$ that you seek.

I did just that and ended up with $\vec r=\langle 1,2,4\rangle^T$ after having reduced the size of the vector. If this is correct, the problem has no solutions since then the line through $P$ never intersects the first line, as it is parallel to it.


text:
\subsection{Monotonic splines}
\label{sec:splines_monotonic}
The spline formulation, presented in Section \ref{sec:splines} to globally enforce a bound constraint of the form $f \ge a$  may be extended easily to enforce monotonicity constraints or other linear inequalities. For example, if $C$ is a first-order (backward or forward) finite difference matrix relating neighboring spline values, then monotonicity is enforced globally by sampling values of the knots $\bm{\xi}$ subject to the constraint
\begin{equation}
C\bm{\xi} \geq \mathbf{0};
\end{equation}
see \citet{lopez2018} or \citet{maatouk}.
This inequality is also used in the rejection sampler of Section \ref{sec:spline_sampling} as a constraint to identify the MAP estimate to increase the sampling efficiency.
Bound and monotonicity constraints can be enforced simultaneously by requiring both $\bm{\xi} \geq \mathbf{b}$ and $C\bm{\xi} \geq \mathbf{0}$ in the sampling, though the acceptance ratio drops substantially with combined constraints.

config: arxiv


TITLE: Question about finding the centre of mass.
QUESTION [0 upvotes]: I'm working through Paul's Online Notes, and I'm on this page.  I don't understand where the equations of moments come from.  Can someone explain to me in the simplest terms how to derive these?
He explains the equations as

the tendency of the region to rotate about the x and y-axis respectively.

So if I wanted to do this with some function $f$, how does the equation
$M_x=\frac{1}{2}\int_a^b f(x)^2 dx$
express the function $f$'s tendency to rotate around the x-axis?

REPLY [3 votes]: Suppose that you have a point mass $m$ located at the point $(x,y)$. By definition its moment about the $x$-axis is $my$, and its moment about the $y$-axis is $mx$. In other words, the moment about an axis is the signed product of the mass and its distance from the axis. Think of a seesaw (teeter-totter): the further you sit from the centre, and the more massive you are, the stronger the tendency for your end to go down.
The nice thing about moments is that they’re additive: the moment of a whole system of point masses about a given axis is the sum of the individual moments. Of course, when you have a whole region of point masses, you can’t simply: you have to integrate. Pick any $x$ between $a$ and $b$ and look at the vertical strip of mass there between $g(x)$ and $f(x)$. A point mass $m$ anywhere along that strip has moment $mx$ about the $y$-axis. The total mass of the strip is $$dm=\rho\big(f(x)-g(x)\big)dx\;,$$ the density times the ‘area’ of the strip, so its moment about the $y$-axis is $$dM_y=xdm=\rho x\big(f(x)-g(x)\big)dx\;.$$ To get the moment about the $y$-axis of the entire region, you have to ‘add up’ (integrate) these individual moments over the entire range of $x$ values, from $a$ to $b$: $$M_y=\rho\int_a^b\big(f(x)-g(x)\big)dx\;.$$
It’s a little harder to see where the formula for $M_x$ comes from, though the basic idea is actually the same. The problem is that each point mass along that vertical strip at $x$ has a different $y$-coordinate, so to get $dM_x$ you’re trying to ‘add up’ $y\;dm$’s with different values of $y$. It turns out that since the density is constant, you can pretend that the entire mass of the strip is concentrated at its midpoint, i.e., at the point $$\left(x,\frac{f(x)+g(x)}2\right)\;,$$ with $y$-coordinate $\bar y=\frac12\big(f(x)+g(x)\big)$. (You’ll be able to show this after you learn about multiple integrals.) The mass of the strip is the same as before, $dm=\rho\big(f(x)-g(x)\big)dx$, so its contribution to the moment about the $x$-axis is $$\begin{align*}
dM_x&=\bar y dm\\
&=\frac12\big(f(x)+g(x)\big)\cdot\rho\big(f(x)-g(x)\big)dx\\
&=\frac12\rho\Big(\big(f(x)\big)^2-\big(g(x)\big)^2\Big)dx\;,
\end{align*}$$
and the actual moment about the $x$-axis is $$M_x=\frac12\rho\int_a^b\Big(\big(f(x)\big)^2-\big(g(x)\big)^2\Big)dx\;,$$ the result of ‘adding up’ the moments for the vertical strips between $x=a$ and $x=b$.


TITLE: A variation of a 'leaky bucket' problem
QUESTION [2 upvotes]: I'm looking at a 'leaky bucket' problem with a few differences, the system is set up as follows:

the tank has a practically unlimited capacity
water is added with a known flow rate
periodically the drain valve is adjusted so that all of the water is drained over say 10 seconds
we are interested to know how much water is contained at time $t$

To give a concrete example, lets assume the water is added with a rate of 10L per second, every 4 seconds the drain is adjusted. When the drain is adjusted it is adjusted to $\frac{1}{10}$ of current volume of water.
Initially I started with the following diff. eq. $dV=I-\frac{1}{10}V$ where $V$ stands for volume and $I$ is the incoming flow rate in L per sec. The issue is that in this equation the drain is constant with respect to volume, whereas in the actual problem the drain is adjusted periodically and in between those events the water is draining with a set flow rate.
I have a feeling it involves using a $\sin$ function to model the period, but I'm not certain as to how to go about it.

REPLY [2 votes]: Say volume of water at $t = 0$ is $V_0$.
Volume of water at time $t = 4n$ where $n \in 0, \mathbb{Z+}$ is then given by
$\displaystyle V_{4n} = \bigg(\frac{3}{5}\bigg)^n V_0 + 10I \bigg(1 - \bigg(\frac{3}{5}\bigg)^n \bigg)$
And volume of water at time $4n+t$ is given by,
$\displaystyle V_{4n+t} = V_{4n} \bigg(1 - \frac{t}{10}\bigg) + I \, t$ where $0 \leq t \leq 4$
Here is how the expression comes (wherever I write $t$, it is $0 \leq t \leq 4$) -
$V_t = V_0 - \frac{V_0}{10}t + I \, t$
$V_4 = V_0 - \frac{4V_0}{10} + 4I = \frac{3}{5}V_0 + 4I$
$\displaystyle V_{4+t} = V_4 - \frac{V_4}{10}t + I \, t$
i.e. $\displaystyle V_{4+t} = \frac{3}{5}V_0 + 4I - \frac{t}{10}(\frac{3}{5}V_0 + 4I) + I \, t$
i.e. $\displaystyle V_{4+t} = \frac{3}{5}V_0 (1 - \frac{t}{10})+ (4 + \frac{3t}{5})I$
$\therefore \displaystyle V_{8} = (\frac{3}{5})^2 V_0+ 4(1 + \frac{3}{5})I$
Similarly, $\displaystyle V_{12} = (\frac{3}{5})^3 V_0 + 4(1 + \frac{3}{5} + (\frac{3}{5})^2)I$
At $t = 4n$, the geometric series $\, 1 + \frac{3}{5} + ... + (\frac{3}{5})^{n-1} = \displaystyle \frac{1 - (\frac{3}{5})^{n}}{1 - \frac{3}{5}} = \frac{5}{2} (1 - (\frac{3}{5})^{n})$
And this leads to expression at time $4n$ and $4n + t$.


TITLE: $\forall x\in\mathbb{R}\enspace \enspace, \forall y\in\mathbb{R} : (x\ge y)$ or $(x\le y)$ proof?
QUESTION [1 upvotes]: Well, it's obvious that it's true.
but our teacher just proved it by saying "for each $x$ and $y$ we will input, we will receive true element".
My question is - How can I prove it in the best way?
I did the contradiction proof and said
$\exists (x,y)\in\mathbb{R} : (x<y)$ AND $(x>y)$.
I input one time $x < y$ , $x=y$ and $x>y$ numbers $($like $2,3$ , $2,2$ , $2,1$$)$
and showed in truth table that its always false which means the main claim is true.
Is it ok to do so or there are better proofs?
Thank you very much

REPLY [2 votes]: It comes from definition of Real numbers - Cauchy sequences, Dedekind cuts or other. In any way aim is to have so called The complete ordered field.

Dedekind cuts:

We take rational numbers $\mathbb{Q}$ and consider its $\alpha$ subsets with following properties:
$$\begin{array}{}
\alpha \neq \emptyset, \alpha \neq \mathbb{Q} \\
q \in \mathbb{Q} \land p \in \alpha \land q<p \Rightarrow q \in \alpha \\
p \in \alpha \Rightarrow \exists r \in \mathbb{Q}, p<r
\end{array}$$
Such subsets are called "cuts". We define for cuts "$\alpha < \beta$",  when $\alpha$ is proper subset of $\beta$. Now can be obtained all desired properties of order, including least upper bound property.

Cauchy sequences:

Again we consider rational numbers $\mathbb{Q}$ with metric $|x-y|$ and look at all Cauchy sequences. We called 2 Cauchy sequence equivalent, when they have same limit. By this we divide $\mathbb{Q}$ into equivalence classes, which we called real numbers. Order for equivalence classes $(x_n) \leqslant (y_n)$ iif $(x_n)$ is equivalent to $(y_n)$ or $\exists N, n>N, x_n \leqslant y_n$. Now again is possible to have all order properties.


TITLE: PDE - 1D Heat Diffusion Problem
QUESTION [2 upvotes]: I am stuck on this problem I have to do for my PDE class.
Consider heat diffusion in a 1D object between x=0 and x=L. It can be shown that the heat flux at the surface is $-k\nabla{}T*\hat{n}$ i.e. its is proportional to the normal derivative at the surface. In our case, we specify adiabatic boundary conditions as follows : $$u_x(0,t)= 0 $$ $$ u_x(L,t)=0$$ If the initial temperature is $f(x)$ solve for the temperature as a function of x and t.
To be honest, I don't really know where to start. So far I have tried to solve the diff. eq : $$  \frac{\partial u}{\partial t} =-k\nabla{}u*\hat{n}$$ using separation of variables but I got nowhere. Could anyone help me understand how to properly start this problem ? Thank you !

REPLY [1 votes]: I will omit the derivation of the heat equation, but I believe you should have
$$
\frac{\partial u}{\partial t} = k\Delta u.
$$
The idea is that for any small region $V$ we have
$$\frac{d}{dt} \int_V u = \int_{\partial V} -k\nabla u \cdot {\bf n}dS.$$
This equation is simply saying that the energy change in the region $V$ (LHS) is given by the flux through the boundary (RHS). Using the divergence theorem on the right hand side you can derive the heat equation.
Now, in 1D, on the interval $[0,L]$ we have simply
$$u_t = k u_{xx}.$$
Assume that $u(x,0)= f(x)$ and $u_x(0,t)=u_x(L,t)=0$. You can use separation of variables from here: assume that $u(x,t) = X(x)T(t)$ and see that
$$XT' = kX''T.$$
Rearranging we have
$$\frac{T'}{kT} = \frac{X''}{X} = -\lambda.$$
This gives us two ODEs:
$$T'+k\lambda T=0$$
and
$$X''+\lambda X=0.$$
Note the boundary values: $u_x(0,t) = X'(0)T(t)=0$ and $u_x(L,t)=X'(L)T(t)=0$.
If $T(t) =0 $ then $u\equiv 0$ and we only recover the trivial solution. So, we may assume that $X'(0)=X'(L)=0$. I believe the rest of the problem should be calculations from ODEs.


TITLE: spaces of probability measures on a Polish space and the convergence
QUESTION [1 upvotes]: I want to read the topic "spaces of probability measures on a Polish space and the convergence". What is the best reference for that ?

REPLY [2 votes]: Converting Ofer Zeitouni's comment to an answer:

Billingsley, Patrick.  Convergence of Probability Measures.
Ethier, Stewart and Kurtz, Thomas.  Markov Processes: Characterization and Convergence.


TITLE: Multiplication of sums
QUESTION [0 upvotes]: In this set of online math notes the following is stipulated:
$$\left(\sum_{n=0}^\infty a_n\right)\left(\sum_{n=0}^\infty b_n\right)  =\sum_{n=0}^\infty (c_n) $$
Where apparently $c_n = \sum_{i=0}^n a_ib_{n-i}$. Replacing $c_n$ into above equation we get
$$\left(\sum_{n=0}^\infty a_n\right)\left(\sum_{n=0}^\infty b_n\right)  =\sum_{n=0}^\infty \left(\sum_{i=0}^n a_ib_{n-i}\right) $$
I tried to expand a few terms of $c_n$ using sympy:
$$ {c}_{0} = \sum_{i=0}^{0} {a}_{i} {b}_{- i} = {a}_{0} {b}_{0} $$
$$ {c}_{1} = \sum_{i=0}^{1} {a}_{i} {b}_{- i + 1} = {a}_{0} {b}_{1} + {a}_{1} {b}_{0} $$
$$ {c}_{2} = \sum_{i=0}^{2} {a}_{i} {b}_{- i + 2} = {a}_{0} {b}_{2} + {a}_{1} {b}_{1} + {a}_{2} {b}_{0} $$
Then I ran this example on a verifiable multiplication of sums as follows:
$$(2+x+2x^2)(3-5x+x^2) = 6-7x+3x^2-9x^3+2x^4$$
By setting
\begin{equation*}
a_0 = 2, a_1 = x, a_2 = 2x^2
\end{equation*}
and
\begin{equation*}
b_0 = 3, b_1 = -5x, b_2 = x^2
\end{equation*}
Following the online reference this should equal
$$\sum_{n=0}^2 c_n =c_0+c_1+c_2 =a_0b_0+a_0b_1+a_1b_0+a_0b_2+a_1b_1+a_2b_0 $$
However when you substitute in the values you get
$$ 3 x^{2} - 7 x + 6 $$
Which is not equal to $6-7x+3x^2-9x^3+2x^4$.
I have a feeling the notes are wrong here. Could anyone confirm this? THanks.

REPLY [2 votes]: The added terms, $-9x^3+2x^4$, will come from $c_3$ and $c_4$, in particular $a_1b_2+a_2b_1$ for $c_3$ and $a_2b_2$ for $c_4$.  You have truncated the $c$ sum too early to see them.

REPLY [1 votes]: The first identity should read $\sum_{n\ge 0}a_n x^n\sum_{n\ge 0}b_n x^n=\sum_{n\ge 0}(\sum_{i= 0}^n a_i b_{n-i})x^n$, or even better $\sum_{j\ge 0}a_j x^j\sum_{k\ge 0}b_k x^k=\sum_{n\ge 0}(\sum_{i= 0}^n a_i b_{n-i})x^n$. The proof is simple: just find all contributions to the $x^n$ coefficient on the left-hand side. When multiplying polynomials, which have finite degrees, all sufficiently late terms in each sequence are $0$. In this case $a_n=0$ for $n>2$, $b_n=0$ for $n>2$, $c_n=0$ for $n>4$.


TITLE: injection with sobolev space
QUESTION [0 upvotes]: Let $\Omega $ be a bounded open subset of $R^n,\; n\ge 1.$
I m asking about the existence of a subregion $\omega\subset \Omega$ such that the map $y\to y|_\omega $ from $H^2(\Omega)$ into $L^\infty(\omega)$ is continuous?.

REPLY [2 votes]: I think it would help to specify what you mean by subregion $\omega \subset \Omega$.  In general, it is absolutely true for $\omega$ a lower dimensional object.  For example, when $n=2$ and $u \in H^1(\Omega)$, you can expect that for most (Lebesgue almost every) slices the restriction $u|_\omega$ is absolutely continuous (and therefore continuous and bounded).
If you want a subregion of full measure, then the dimension and exponent are important.  A nice example is a function like $|x|^a$, whose derivative is $a|x|^{a-1}$ and second derivative is $a(a-1)|x|^{a-2}$.  Computing the $L^2$ norm of this second derivative near zero, we find that we require
$\int_0^1 r^{2(a-2)}\;r^{n-1}dr <\infty$
which is equivalent to saying $2(a-2)+n-1>-1$, so that we should have $a>2-\frac{n}{2}$.  Therefore, if $n > 4$, the function $|x|^{-\epsilon}$ is in $H^2(B)$ but not in $L^\infty(B)$.  Then, you can define
$u(x):= \sum_n \frac{1}{2^n} |x-x_n|^{-\epsilon}$,
where $(x_n)_n$ is a dense set and obtain a function which is unbounded in every open set and still in $H^2$, by the above calculation.
Also, YangMills comment is right, since we see when $n=4$ it should still be possible, but one should use logarithms.


text:
\begin{document}

\maketitle

\begin{abstract}
In this paper, we consider an online enrichment procedure using the Generalized Multiscale Finite Element Method (GMsFEM) in the context of a two-phase flow model in heterogeneous porous media. The coefficient of the elliptic equation is referred to as the permeability and is the main source of heterogeneity within the model. The elliptic pressure equation is solved using online GMsFEM, and is coupled with a hyperbolic transport equation where local conservation of mass is necessary. To satisfy the conservation property, we aim at constructing conservative fluxes within the space of multiscale basis functions through the use of a postprocessing technique. In order to improve the accuracy of the pressure and velocity solutions in the online GMsFEM we apply a systematic online enrichment procedure. The increase in pressure accuracy due to the online construction is inherited by the conservative flux fields and the desired saturation solutions from the coupled transport equation. Despite the fact that the coefficient of the pressure equation is dependent on the saturation which may vary in time, we may construct an approximation space using the initial coefficient where no further basis updates follow. Numerical results corresponding to four different types of heterogeneous permeability coefficients are exhibited to test the proposed methodology.
\end{abstract}

\keywords{online enrichment  \and flows in heterogeneous media \and Generalized Multiscale Finite Element Method\and Postprocessing}

\section{Introduction}
A large number of problems of fundamental and practical significance are described by partial differential equations with coefficients that vary over a wide range of length scales. For example, composite materials, porous media and turbulent transport in high Reynolds number flows are models of this type. The heterogeneity and high-contrast properties of the coefficients cause significant difficulty in analyzing these types problems. In this paper, we consider a two-phase flow model in which the so-called permeability coefficient is assumed to be highly heterogeneous. Solving this type of model problem on a fine scale that sufficiently captures the underlying behavior of the heterogeneity may become prohibitively expensive. As a result, methods that aim toward effectively reducing the dimension of the associated fine-scale system(s) have been a topic of continued interest in recent decades. For example, upscaling procedures (see, e.g., \cite{chen2003coupled,durlofsky1991numerical,wu2002analysis}) and multiscale methods (see, e.g., \cite{efendiev2011multiscale,jenny2003multi,wheeler2012multiscale,efendiev2009multiscale,hou1997multiscale} ) are approaches that have been shown to offer effective alternatives to direct fine-scale computations.
For upscaling, one derives a set of localized problems in which averaged quantities may be maintained while solving a lower dimensional global problem on a coarse grid. However, this type of approach may diminish important fine-scale information that has strong effect on the solution behavior. Multiscale methods, on the other hand, hinge on the the independent construction of a set of multiscale basis functions that are used to span a coarse-grid solution space. The coarse-grid discretization parameter may be much larger than the characteristic scale of heterogeneous coefficient, however, the multiscale basis functions inherently include the fine-scale information of the underlying heterogeneity of the medium.

In order to model multi-phase flow, local mass conservation for the fluid velocity fields is required. This requirement has motivated a variety of mass-conservative approaches, such as multiscale finite volume methods \cite{cortinovis2014iterative,jenny2003multi,lunati2004multi}, mixed multiscale finite element methods \cite{aarnes2004use,aarnes2008mixed,chen2003mixed,chung2015mixed}, mortar multiscale methods \cite{arbogast2007multiscale, peszynska2005mortar,peszynska2002mortar}, discontinous Galerkin (DG) methods \cite{du2018adaptive,kim2013staggered,cockburn2002local}, and postprocessing methods \cite{odsaeter2017postprocessing,bush2013application}. In this paper, we use a global continuous Galerkin (CG) method. An advantage of CG multiscale formulation is the relative ease of implementation. On the other hand, a CG solution does not automatically satisfy local conservation, which is essential in our model problem. In order to address this limitation, we adpot an analogous postprocessing technique from \cite{bush2014application}. In particular, after obtaining a multiscale solution, we solve an independent set of local auxiliary problems in order to obtain the locally conservative fluxes.

In terms of the standard Multiscale Finite Element Method (MsFEM), there are two main shortcomings. The first one is that only one basis in each local neighborhood may not be sufficient to guarantee an accurate approximation, especially when there are long channels and non-separable scales in the permeability field. The second limitation is that we often assume that the local boundary conditions are linear along the edges of coarse blocks, which may create a mismatch between multiscale solution and fine-scale solution on the coarse block boundaries. One such technique that may be used to reduce the effect of boundary terms is oversampling \cite{efendiev2000convergence,hou1997multiscale}. Oversampling involves the enlargement of the local computing regions in order to address the linear boundary values. A more recent method that serves to improve the accuracy of MsFEM is the Generalized Multiscale Finite Element Method (GMsFEM) \cite{efendiev2013generalized}.
GMsFEM is a flexible general framework that generalizes MsFEM by systematically enriching the coarse spaces. In particular, more basis functions are added to the initial approximation space in order to improve the accuracy of the multiscale solution. The creation of GMsFEM solution spaces often involves the construction of snapshot, offline, and online spaces \cite{chung2015residual,chung2017online} in order to streamline the procedure for repeated basis function computations. In order to construct an offline multiscale space, some well-designed local spectral problems are solved in order to obtain a set of basis functions that are independent of global information such as source terms and boundary conditions. These local problems are motivated by the convergence analysis, which offers a convergence rate of $1/\Lambda$, where $\Lambda$ is the smallest eigenvalue whose modes are excluded in the multiscale space. If we increase the number of offline bases to a certain number, the error decay will diminish, and in \cite{CHUNG201669}, it is shown that a good approximation from the reduced model can be expected only if the offline information is a good representation of the problem.  Consequently, an online enrichment procedure is essential if the offline bases are not sufficiently accurate. The main idea in this paper is to enrich the offline space by incorporating a new set of basis functions in order to obtain a significant error decay. In \cite{chan2016adaptive,chung2015residual}, the authors propose an online construction resulting from the associated offline space. In consideration of fact that the offline bases are obtained independently through a set of local problems, one may seek to construct a set of bases that contain some global information. Based on this idea, we use residual-driven basis functions which are computed through a set of local problems. The analysis in \cite{chung2017online,chung2015residual} shows that the error decay is proportional to $1-\Lambda$.

The rest of paper is organized as follows. In Section 2 we introduce the model problem and the corresponding solution algorithm. In Section 3, we describe the Generalized Multiscale Finte Element Method (GMsFEM) and the construction of the online solution space. The post-processing technique that is used to ensure the local conservation of mass property is reviewed in Section 4. In Section 5 we offer a variety of numerical results to illustrate the effectiveness of the proposed methodology.

\section{Model problem}

\subsection{Two-phase model}
In this paper, we consider the dynamics of  the movements of two immiscible fluids in a heterogeneous oil reservoir constrained in a domain $\Omega$. In particular, we model scenario where water is discharged to replace trapped oil in a saturated subsurface. Under the assumptions that the environment is gravity-free, capillary pressure is not included, and that two fluids fill the pore space we can apply Darcy's law combined with a statement of conservation of mass. The principle equations of the flow may then be stated as follows:
\begin{eqnarray}
\nabla\cdot \boldsymbol{v}=q, ~~ \text{where} ~~ \boldsymbol{v}=-\lambda(S)k(x)\nabla p\label{elliptic}
\end{eqnarray}
\begin{eqnarray}
\frac{\partial S}{\partial t}+\nabla \cdot(f(S) \boldsymbol{v})=q_{w}, \label{transport}
\end{eqnarray}
where $p$ is the pressure, $\mathbf{v}$ is the Darcy velocity, $S$ is the water saturation, $q,q_w$ are any external forces and $k(x)$ is the heterogeneous permeability coefficient. The total mobility $\lambda(S)$ and the flux function $f(S)$ are respectively given by:
\begin{eqnarray*}
\lambda(S)=\frac{k_{r w}(S)}{\mu_{w}}+\frac{k_{r o}(S)}{u_{e}}, \quad f(S)=\frac{k_{r w}(S) / \mu_{w}}{\lambda(S)}
\end{eqnarray*}
where $k_{r,j},j=w,o$, is the relative permeability of the phase $j$.
\subsection{Solution algorithm}
In Table \ref{two-phase algorithm}, we display the algorithm that is used to solve the two-phase model in Eqs. \eqref{elliptic} and \eqref{transport}. In order to solve for the unknown saturation $S$, we first split the time interval into a set of specified subintervals. $S$ is initialized by $S_0$ and then solved by a series of iterations that are included in Table \ref{two-phase algorithm}. More specifically, we use $S_{n-1}$ in \eqref{elliptic} to obtain $p_n$ and $v_n$. Then we solve \eqref{transport} using the new flux $\mathbf{v}_n$ to obtain $S_n$.\\
\begin{table}[!htbp]
    \begin{tabular}{c l}
    \hline
      &Two-phase algorithm  \\
       \hline
    \textbf{Input} & $S_{n-1}$ obtained in previous time step\\
      \textbf{Output}& $S_{n}$ \\
      &1. Solving \ref{elliptic} to get $p_n$ and $\mathbf{v}_n$\\
      &2. Using $\boldsymbol{v}_n$ and $S_{n-1}$  in \ref{transport} to get $S_{n}$ \\
      \hline
      \caption{Two-phase algorithm}
      \label{two-phase algorithm}
      \end{tabular}
\end{table}

To solve \eqref{transport}, we integrate over the time interval $[t_{n-1},t_n]$ and a control volume  $C_z\subset \Omega$ to obtain
\begin{eqnarray}
\text{meas}(C_z)(S_{z,n}-S_{z,n-1})+\Delta t\int_{\partial_{C_z}}\mathbf{v}\cdot \mathbf{n} f(S_{z,n-1}) ~ dl=\Delta t\int_{C_z} q_w ~ dx,
\end{eqnarray}
where we have neglected the error terms, and we use
\begin{eqnarray}
S_{z,n}\approx\dfrac{1}{meas(C_z)}\int_{C_z}S(x,t_n) ~ dx.
\end{eqnarray}
We use $meas(A)=\int_{\Omega} 1_{A}d x$ with $1_A=1$ when $x\in A$ while 0 elsewhere.
To evaluate the term $\int_{\partial_{C_z}} \mathbf{v} \cdot \mathbf{n} f(S_{z,n-1}) \, dl$, we use an upwinding scheme. A review of upwinding on a rectangular mesh can be in \cite{thomas2013numerical}, for example. It is imperative that the numerical approximation of $\mathbf{v}$ satisfies the following local conservation property. In particular, it is desirable to have
\begin{eqnarray}
\int_{\partial{C_z}}\mathbf{v}\cdot \mathbf{n} ~d l=\int_{C_z} q ~d x.
\end{eqnarray}
There are two main ways to obtain the desired quantities $(\mathbf{v},p)$. The first one is to simultaneously solve the first order system \eqref{elliptic}. For example, one may apply the mixed finite element formulation \cite{chan2016adaptive}. In this paper, we consider the alternative of transforming \eqref{elliptic} into a second order equation that  governs the pressure  $p$. The approximation of $\mathbf{v}$ is calculated using the relation $\mathbf{v}=-\lambda(S)k(x)\nabla p$, and a postprocessing procedure follows for local conservation. Since it is computationally expensive to apply the postprocessing procedure on the fine-scale solution, we instead use the Generalized Multiscale Finite Element Method (GMsFEM), which will be introduced in the next section.

\section{Generalized multiscale finite element method}
\subsection{Preliminaries}
We fix our attention to the following second order elliptic problem
\begin{eqnarray}
\begin{aligned}
-\operatorname{div}(\lambda k(x) \nabla p)&=q \quad \text { in } \Omega \\
p&=p_{D} \quad \text { on } \Gamma_{D} \\
-\lambda k \nabla p \cdot \mathbf{n}&=g_{N} \quad \text { on } \Gamma_{N}\label{model}
\end{aligned}
\end{eqnarray}
where $k(x)$ is a highly heterogeneous field with high contrast. In practice, we assume that there is a positive  constant $k_{min}$ such that $k(x)\geq k_{min}\geq0$, while $k(x)$ can vary widely (i.e.,
$k_{max}/k_{min}$ is very large, for example $10^5$). Four examples of $k(x)$ that are considered in this paper are offered in Figure \ref{medium}. All permeability fields in the figure are plotted on the log scale. Additionally, $\lambda$ is a known mobility coefficient, $q$ denotes any external forcing, and $p$ is an unknown pressure field satisfying Dirichlet and Neumann boundary conditions given by $p_D$ and $g_N$, respectively. Here $\Omega$ is a convex polygonal and two dimensional domain with boundary $\partial{\Omega}=\Gamma_{D}\cup \Gamma_{N}$.

We consider a function in $H^1(\Omega)$ whose trace on $\Gamma_{D}$ coincides with the given value $p_D$; we
denote this function also by $p_D$. The variational formulation of \eqref{model} is stated as follows. We find $p\in H^1({\Omega})$ with $(p-p_D)\in H^1_{\Omega}=\{w\in H^1_0(\Omega):w|_{\Gamma_{D}}=0 \}$ such that
\begin{eqnarray}
a(p,v)=F(v)-\left\langle g_N,v\right\rangle_{\Gamma_{N}} ~~\text{for all } ~~v\in H^1_D \label{varational}
\end{eqnarray}
where
\begin{eqnarray*}
\begin{aligned}
a(p, v)&=\int_{\Omega} \lambda k(x) \nabla p(x) \nabla v(x)~ d x,\\
F(v)&=\int_{\Omega}q(x)v(x) ~d x, ~~\text{and}\\
\left\langle g_{N}, v\right\rangle_{\Gamma_{N}}&=\int_{\Gamma_{N}} g_{N}(x) v(x)~ d l.
\end{aligned}
\end{eqnarray*}
\begin{figure}[!htbp]
\centering
\subfigure[$\kappa_1(x)$]
{ \includegraphics[width=0.45\textwidth]{graphs/s/k1.png}}
\subfigure[$\kappa_2(x)$]
{ \includegraphics[width=0.45\textwidth]{graphs/s/k2.png}}
\subfigure[$\kappa_3(x)$]
{ \includegraphics[width=0.45\textwidth]{graphs/s/k3.png}}
\subfigure[$\kappa_4(x)$]
{ \includegraphics[width=0.45\textwidth]{graphs/s/k4.png}}
\caption{Examples of heterogeneous permeability fields; all plots are on the log scale}\label{medium}
\end{figure}
 \begin{figure}[ht]
    \centering
 {\includegraphics[width=3in]{graphs/coarse_neighborhood.png}}
 \caption{Discretization of $\Omega$ into $\mathcal{T}_h=\cup \tau$. Here $\omega_z=\cup_{i=1}^4\tau_i$ is the supp($\chi_z$).}\label{figure:coarse}
 \end{figure}

In order to implement a finite element approximation of \eqref{varational}, we
let $\mathcal{T}^{h}$ denote a partition of the domain $\Omega$ into fine elements. Here, $h>0$ is used to denote the fine-grid mesh size.
 The coarse partition, $\mathcal{T}^{H}$ of the domain $\Omega$, is formed such that each element in $\mathcal{T}^{H}$ is a connected union of fine-grid blocks. More precisely, $\forall K_{j} \in \mathcal{T}^{H}$, $ K_{j}=\bigcup_{F\in I_{j} }F$ for some $I_{j}\subset \mathcal{T}^{h}$. The quantity $H>0$ is the coarse mesh size. In this paper we consider the case of rectangular coarse elements, yet the methodology can be used with general coarse elements. An illustration of the mesh notations is shown in the Figure \ref{figure:coarse} (\textbf{the notation in the illustration does not match the notation used below. For example $\omega_z$ is used in the figure for a neighborhood, whereas $D_i$ is used below for the neighborhood}). We denote the interior nodes of $\mathcal{T}^{H}$ by $x_i, ~~i=1,\cdots,N_{\text{in}}$,
 where $N_\text{in}$ is the number of interior nodes. The coarse elements
 of $\mathcal{T}^{H}$ are denoted by $K_j, ~~j=1,2,\cdots,N_e$, where $N_e$ is the number of coarse elements. We define the coarse neighborhood of the nodes $x_i$ by $D_i:=\cup\{K_j\in T^{H}:x_i\in \overline{K_j}\}$.

\subsection{GMsFEM for pressure equation}
 In this paper, we will apply the GMsFEM to solve nonlinear parabolic equations. The
 method is motivated by the finite element framework. First, a variational formulation is defined. Then we construct some multiscale basis functions.
  Once the fine grids are given, we can compute the fine-grid solution. Let $\gamma_1,\cdots,\gamma_n$ be the standard finite element basis, and define $V_f=\text{span}\{\gamma_1,\cdots,\gamma_n\}$ to be the
 fine space.
 We obtained the fine solution denoted by $p_h$ by solving
 \begin{eqnarray}
a(p_h,v_h)=F(v_h)-\left\langle g_N,v_h\right\rangle_{\Gamma_{N}} \text{for all }v_h\in V_f \label{fine}
 \end{eqnarray}
  The construction of multiscale basis functions follows two general steps. First, we construct snapshot basis functions in order to build a set of possible modes of the solutions. In the second step, we construct multiscale basis functions with a suitable spectral problem defined in the snapshot space. We take the first few dominated eigenfunctions as basis functions. Using the multiscale basis functions, we obtain a reduced model.
More specifically, once the coarse and fine grids are given, one may construct the multiscale basis functions to approximate the solution
of (\ref{varational}). To obtain the multiscale basis functions, we first define the snapshot space. For each coarse neighborhood $ D_{i}$, define $J_h( D_{i})$ as the set of the fine nodes of $T^{h}$ lying on $\partial D_{i}$ and denote
the its cardinality by $L_i \in \mathbb{N}^{+}$. For each fine-grid node $x_j \in J_h( D_{i})$, we define a fine-grid function $\delta_{j}^{h}$ on $J_h( D_{i})$ as $\delta_{j}^{h}(x_k)=\delta_{j,k}$. Here
 $\delta_{j,k}=1$ if $j=k$ and $\delta_{j,k}=0$ if $j\neq k$. For each $j=1,\cdots, L_i$, we define the snapshot basis functions $\psi_{j}^{(i)}$ ($j=1,\cdots,L_i$) as the solution of the following system\\
\begin{eqnarray}
\begin{aligned}
   -\nabla\cdot\left(\kappa \nabla \psi_{j}^{(i)}\right) &=0 \quad \text { in } D_{i} \\
   \psi_{j}^{(i)} &=\delta_{j}^{h} \quad \text { on } \partial D_{i}.\label{snap_basis}
\end{aligned}
\end{eqnarray}
 The local snapshot space $V_{\text { snap }}^{(i)}$ corresponding to the coarse neighborhood $ D_{i}$  is defined as follows
  $V_{snap}^{(i)}:=$ \text{span}$\{\psi_{j}^{(i)}:j=1,\cdots,L_{i}\}$ and the snapshot space reads $V_{\text {snap}} :=\bigoplus_{i=1}^{N_{\text {in}}} V_{\text {snap}}^{(i)}$, where $N_{\text {in}}$ is the total number of coarse neighborhood.

  In the second step, a dimension reduction is performed on $V_{\text {snap}}$.
  For each $i=1,\cdots, N_{\text {in}}$, we solve the following spectral problem:
  \begin{eqnarray}
      \int_{D_{i}} \kappa \nabla \phi_{j}^{(i)} \cdot \nabla v=\lambda_{j}^{(i)} \int_{D_{i}} \hat{\kappa} \phi_{j}^{(i)} v \quad \forall v \in V_{\text {snap}}^{(i)}, \quad j=1, \ldots, L_{i}\label{eigen}
  \end{eqnarray}
  where $\hat{\kappa} :=\kappa \sum_{i=1}^{N_{i n}} H^{2}\left|\nabla \chi_{i}\right|^{2}$ and $\{\chi_{i}\}_{i=1}^{N_{i n}}$ is a set of partition of unity that solves the following system:
  \begin{eqnarray*}
  \begin{array}
  {rlrl}{-\nabla \cdot\left(\kappa \nabla \chi_{i}\right)} & {=0} & {} & {\text { in } K \subset D_{i}} \\ {\chi_{i}} & {=p_{i}} & {} & {\text { on each } \partial K \text { with } K \subset D_{i}} \\ {\chi_{i}} & {=0} & {} & {\text { on } \partial D_{i}}
  \end{array}
  \end{eqnarray*}
  where $p_i$ is some polynomial functions and we can choose linear functions for simplicity.
  Assume that the eigenvalues obtained from (\ref{eigen}) are arranged in ascending order and we may use the first $1<l_i \leq L_{i}$ (with $l_{i} \in
  \mathbb{N}^{+}$) eigenfunctions (related to the smallest $l_i$ eigenvalues) to
  form the local multiscale space $V_{\text{off}}^{(i)}:=$ snap$\{\chi_{i}\phi_{j}^{(i)}:j=1,\cdots,L_{i}\}$. The mulitiscale space $V_{\text{off}}^{(i)}$ is the direct sum of the local mulitiscale spaces, namely
  $V_{\text {off}} :=\bigoplus_{i=1}^{N_{\text {in}}} V_{\text {off}}^{(i)}$.
  Once the multiscale space $V_{\text {off}}$ is constructed, we can find the
  GMsFEM solution $p_H$ by solving the following equation
 \begin{eqnarray}
a(p_H,v_H)=F(v_H)-(g_N,v_H)_{\Gamma_{N}} \text{for all }v_H\in V_{\text {off}} \label{coarse}
 \end{eqnarray}
 In the numerical examples, we use $L_z$ to denote $L_i$ for $1\leq i\leq N_{\text {in}}$ since we use same $L_i$ for each $i$.
\subsection{Online enrichment}
We will present the constructions of online basis functions \cite{chung2015residual} in this section. \\
After obtaining the multiscale space $V_{\text {off}}$, one may add some online basis functions based on local residuals.\\
Let $p_H \in V_{\text {off}}$ be the solution obtained in (\ref{coarse}). Given a coarse neighborhood $D_i$, we define
$V_i:=H_0^1(D_i)\cap V_{\text {snap}}$ equipped with the norm
$\|v\|_{V_i}^{2}:=\int_{D_i}\kappa |\nabla {v}|^2$. We also define the local residual operator $R_i: V_i\rightarrow \mathbb{R}$ by
\begin{eqnarray}
    \mathcal{R}_{i}\left(v ; p_H\right) :=a(p_H,v)-F(v)+(g_N,v)_{\Gamma_{N}} \label{loc res}
\end{eqnarray}
 The norm of operator $R_i$, denoted by $\|R_i\|_{V_{i}^{*}}$, gives a measure of the quantity of residual.

  Suppose one needs to add one new online basis $\phi$ into the space $V_i$. The analysis in \cite{chung2015residual} suggests that the required online basis $\phi\in V_i$ is the solution to the following equation
  \begin{eqnarray}
      \mathcal{A}(\phi, v)=\mathcal{R}_{i}\left(v ; p_H^{\tau}\right) \quad \forall v \in V_{i}.\label{online}
  \end{eqnarray}
  We refer to $\tau \in \mathbb{N}$  as the level of the enrichment
  and denote the solution of (\ref{coarse}) by $p_H^{\tau}$.
  Remark that $V_{\text {off}}^{0}:=V_{\text{off}}$. Let $\mathcal{I} \subset\left\{1,2, \ldots, N_{i n}\right\}$ be the index set over some non-lapping coarse neighborhoods. For each $i\in \mathcal{I}$, we obtain a online basis $\phi_i\in V_i$ by solving (\ref{online}) and define
$V_{\text {off}}^{\tau+1}=V_{\text {off}}^{\tau} \oplus \operatorname{span}\left\{\phi_{i} : i \in \mathcal{I}\right\}$.
After that, solve (\ref{coarse}) in $V_{\text {off}}^{\tau+1}$.

\section{Postprocessing GMsFEM solution}\label{postprocess}
In order to obtain the GMsFEM with local conservation property, we apply postprocessing technique after obtain $p_H$. The technique was introduced in \cite{bush2014application}. In this section, a review is presented.\\
This approach is composed of two main steps. The first step is solving an auxiliary boundary value problem element by element.  The next step is called downscaling procedure, which is solving similar boundary value problem in each control volumn using the auxiliary solutions obtained in first step. Derivation of local conservation of mass is presented in \ref{conservation}.
\subsection{Constructing a locally conservative flux}
In particular, we obtain a auxiliary solution denoted by $\tilde{p}_{\tau}$,with
\begin{eqnarray}
    \left\{\begin{array}{ll}
-\nabla \cdot\left(\lambda \kappa(x) \nabla \tilde{p}_{\tau}\right)=q & \text { in } \tau \\
-\lambda \kappa(x) \nabla \tilde{p}_{\tau} \cdot n=\tilde{g}_{\tau} & \text { on } \partial \tau
\end{array}\right.\label{auxiliary}
\end{eqnarray}
 \begin{figure}[!htbp]
    \centering
    \subfigure
       { \includegraphics[width=0.45\textwidth]{graphs/control_volumn.png}}
        \subfigure
       { \includegraphics[width=0.3\textwidth]{graphs/element.png}}
    \caption{Left: $C_z$ is the control volumn associated with the vertex z,where $\partial C_z=E_{z\eta}\cup E_{Z\omega}\cup E_{z\xi}\cup E_{z\gamma}$.
    Right: a finite element $\tau$ is divided into four quadrilaterals $t_z,t_{\omega},t_x,t_{\gamma}$.}\label{control volumn}
\end{figure}
Here, we designate $\partial{\tau}=\cup_{\xi\in v(\tau)} E_{\xi}^{\tau}$, where $E_{\xi}^{\tau}=\partial \tau \cap \partial t_{\xi}$ (i.e. half of each element edge containing the vertex $\xi$.) and $v(\tau)$ is the collection of four vertexes of $\tau$.
Furthermore, we set $\tilde{g}_{\tau}$ as piecewise function on $\partial\tau$ such that
\begin{eqnarray*}
    \int_{E_{\xi}^{\tau}} \tilde{g}_{\tau} \mathrm{d} l=F_{\xi, 1}-Q_{\xi, 1}, \quad \text { for } \xi \in v(\tau)
\end{eqnarray*}
where
\begin{eqnarray}
Q_{\xi, 1}=\int_{\tau} \lambda \kappa \nabla p_{H} \cdot \nabla \Phi_{\xi, 1} \mathrm{d} x \quad \text { and } \quad F_{\xi, 1}=\int_{\tau} q \Phi_{\xi, 1} \mathrm{d} x \label{Q and F}
\end{eqnarray}
The existence and uniqueness of the above problem is stated in \cite{bush2014application}.
\ref{auxiliary} implies
\begin{equation*}
    -\int_{\partial \tau} \lambda \kappa \nabla \tilde{p}_{\tau} \cdot \boldsymbol{n} \mathrm{d} l=-\int_{\partial \tau} \lambda \kappa \nabla p \cdot \boldsymbol{n} \mathrm{d} l,
\end{equation*}
which shows that the solution of \ref{auxiliary} recovers the flux of p (i.e. the true pressure solution) averaged over $\partial \tau$, a local conservation property in each element. We use \ref{auxiliary} as a governing principle to derive the processing technique for calculating a locally conservative flux in each control volumn from $p_H$.
 The elemental calculation is based on discretization of $\tau$ into quadrilaterals $t_{\xi}$,i.e., $\tau=\cup_{\xi\in v(\tau)}t_{\xi}$, each of which yields $t_{\xi}=C_{\xi}\cap\tau$,see the right plot of Figure \ref{control volumn}. We set the local solution space as $\mathcal{V}(\tau)=\text{span}\{\Phi_{\xi,1}\}_{\xi\in v(\tau)}$, where $\Phi_{\xi,1}$ is the multiscale basis function corresponding to the vertex $\xi$ . The numerical solution associated with
 \ref{auxiliary} is to find $\tilde{p}_{\tau,h}\in \mathcal{V}(\tau)$ satisfying
 \begin{equation}
     -\int_{\partial t_{\zeta}} \lambda \kappa \nabla \tilde{p}_{\tau, h} \cdot \boldsymbol{n} \mathrm{d} l=\int_{t_{\zeta}} q \mathrm{d} x, \quad \text { for all } \zeta \in v(\tau)\label{part}
 \end{equation}
 The following four equations result from \ref{part}:
 \begin{eqnarray}
     \begin{aligned}
&q_{z \omega}^{\tau}+q_{z\gamma}^{\tau}=Q_{z, 1}-F_{z, 1}+\int_{t_{z}} q \mathrm{d} x\\
&q_{z\gamma}^{\tau}+q_{x \gamma}^{\tau}=Q_{\gamma, 1}-F_{\gamma, 1}+\int_{t_{\gamma}} q \mathrm{d} x\\
&q_{x \omega}^{\tau}+q_{z \omega}^{\tau}=Q_{\omega, 1}-F_{\omega, 1}+\int_{t_{\omega}} q \mathrm{d} x\\
&q_{x \gamma}^{\tau}+q_{x \omega}^{\tau}=Q_{x, 1}-F_{x, 1}+\int_{t_{x}} q \mathrm{d} x \label{system}
\end{aligned}
 \end{eqnarray}
 where
 \begin{eqnarray}
     \begin{aligned}
q_{x\omega}^{\tau} &=-\int_{E_{x\omega}^{\tau}} \lambda \kappa \nabla \tilde{p}_{\tau, h} \cdot \boldsymbol{n} \mathrm{d} l, & & q_{z \gamma}^{\tau}=-\int_{E_{z \gamma}^{\tau}} \lambda \kappa \nabla \tilde{p}_{\tau,h} \cdot \boldsymbol{n} \mathrm{d} l \\
q_{z\omega}^{\tau} &=-\int_{E_{z\omega}^{\tau}} \lambda \kappa \nabla \tilde{p}_{\tau, h} \cdot \boldsymbol{n} \mathrm{d} l, & & q_{x \gamma}^{\tau}=-\int_{E_{x \gamma}^{\tau}} \lambda \kappa \nabla \tilde{p}_{\tau, h} \cdot \boldsymbol{n} \mathrm{d} l
\end{aligned}
 \end{eqnarray}
 and $E_{\xi\eta}^{\tau}=\partial t_{\xi}\cap \partial t_{\eta}$, for $\xi,\eta=\omega,x,\gamma,z$ and $\xi\neq\eta$.
 Since we actually use linear combination of basis to solve solution, in particular,
 $\tilde{p}_{\tau,h}=\Sigma_{\xi\in v(\tau)} u_{\xi}\Phi_{xi}$ with unknown coefficients $u_{\xi}$,
 \ref{system} can be written in the form of
 $\tilde{A}\tilde{u}=\tilde{f}$ where

 \begin{equation*}
 \tilde{\boldsymbol{A}}_{\zeta \eta}=-\int_{E_{\mathrm{fp}}} \lambda \kappa \nabla \Phi_{\eta,1}+\boldsymbol{n} \mathrm{d} l \quad \text { and } \quad \tilde{\mathrm{f}}_{\zeta}=\int_{t_{\mathrm{\xi}}} q \mathrm{d} x-\int_{E_{\xi}^{\tau}} \tilde{g}_{\tau} \mathrm{d} l
 \end{equation*}
One should note that when $\tau$ is adjacent to $\Gamma_N$, $g_N$ should be taken in account in computing $\tilde{g}_{\tau}$.\\
Since the system actually has smaller dimension than 4, we may add a constant to one entry in $\tilde{A}$ to remove the singularity. The fact that $u$ is not unique is irrelevant since the desired solution is flux as governed by $q_{\xi\eta}^{\tau}$ which is unique.\\
\ref{system} implies that $\tilde{v}_h$ derived from $\tilde{p}_h$ satisfy the desired local conservation property.
\subsection{Downscale procedure}
After the postprocessing in the section 4.1, we have
\begin{equation*}
    \int_{\partial C_{\tau}} \tilde{v}_{h} \cdot \boldsymbol{n} \mathrm{d} l=\int_{\mathcal{C}_{z}} q \mathrm{d} \boldsymbol{x}, \quad \text { for all } C_{z}
\end{equation*}
which can be thought of a statement of compatibility condition in $C_z$.
We can proceed with formulating a boundary problem as follows,
\begin{eqnarray}
    \left\{\begin{array}{l}
-\nabla \cdot\left(\lambda \kappa(\boldsymbol{x}) \nabla \tilde{p}_{G}\right)=q \quad \text { in } C_{z} \\
-\lambda \kappa(\boldsymbol{x}) \nabla \widetilde{p}_{G} \cdot \boldsymbol{n}=\tilde{v}_{h} \cdot \boldsymbol{n} \text { on } \partial C_{z}
\end{array}\right.\label{downscale}
\end{eqnarray}
 Here $\tilde{v}_{h}=\Sigma_{\tau,\tau\cap C_z\neq \emptyset}-\lambda \kappa(x)\nabla\tilde{p}_{\tau,h}$ that is evaluated pointwise on segments of $\partial C_z$ that belongs to $\tau$.\\
  For example, for control volume $C_z$ corresponding to vertex $z$, we obtain $\tilde{v}_{h}$ as follows. One may refer to the left in figure \ref{control volumn}.
  \begin{eqnarray}
      \begin{aligned}
\int_{\partial C_{z}} \tilde{v}_{h} \cdot n d l &=\int_{E_{z\eta}} \tilde{v}_{h} \cdot n d l+\int_{E_{z\omega}} \tilde{v}_{h} \cdot n d l+\int_{E_{z\xi}} \tilde{v}_{h} \cdot n d l+\int_{E_{z\gamma}} \tilde{v}_{h} \cdot n d l \\
&=\left(q_{z\eta}^{\tau_{1}}+q_{z\eta}^{\tau_{2}}\right)+\left(q_{z \omega}^{\tau_{1}}+q_{z\omega}^{\tau_{3}}\right)+\left(q_{z\xi}^{\tau_{3}}+q_{z\xi}^{\tau_{4}}\right)+\left(q_{z \gamma}^{\tau_{4}}+q_{z \gamma}^{\tau_{2}}\right)\\
&=\left(q_{z\eta}^{\tau_{1}}+q_{z\omega}^{\tau_{1}}\right)+\left(q_{z\eta}^{\tau_{2}}+q_{z\gamma}^{\tau_{2}}\right)+\left(q_{z\omega}^{\tau_{3}}+q_{z_{z\xi}}^{\tau_{3}}\right)+\left(q_{z\xi}^{\tau_{4}}+q_{z\gamma}^{\tau_{4}}\right)\\
&=\sum_{j=1}^{4}\left(Q_{z, 1, j}-F_{z, 1, j}\right)+\sum_{j=1}^{4} \int_{t_{z, j}} q \mathrm{d} \mathbf{x}\\
&=\int_{ C_{z}} q \mathrm{d} \mathbf{x}. \label{conservation}
\end{aligned}
 \end{eqnarray}
 where $Q_{z, 1, j},F_{z, 1, j}$ are integrals (refer \ref{Q and F}) in domain $\tau_j$ for corresponding $j$. $\sum_{j=1}^{4}\left(Q_{z, 1, j}-F_{z, 1, j}\right)=0$ is derived by \ref{coarse}.\\
This calculation actually proves the local conservation of $\tilde{v}_{h}$. So the satisfies compatibility condition of \ref{downscale} guarantees the existence of the corresponding solution. Similarly, since our interest only lies in $-\lambda \kappa(x)\nabla\tilde{p}_{C_z}$ in \ref{downscale}, the nonuniqueness of the solution is of no concern.\\

\section{Numerical results}
In this section, we consider four kinds of permeability coefficients which are represented in figure \ref{medium}. $\kappa_1$ and $\kappa_4$ are extracted from the tenth SPE comparative solution project (SPE10), which is commonly used as benchmark permeability field to assess upscaling and multiscale methods. The most distinguishable characteristic of the model is that some layers are highly heterogeneous and contains long channels. Here, $\kappa_1$ is the last layer of the SPE10 dataset while $\kappa_4$ comes from the 36-th layer. It is evident that $\kappa_1$ represents high heterogeneity and both two contains some visible channels. In terms of $\kappa_2$, it is deterministic, high-contrast coefficient with abrupt transitions between regions of low and high permeability. For $\kappa_3$, it comes from fractured porous media, which is characterized by complex fracture distribution and high contrast. Consequently, four examples of permeability exhibit high-contrast features, which can make solving (\ref{model}) a demanding task.

Since the construction of multiscale space is based on the single-phase flow, i.e. choosing $\lambda(S)=1$ in (\ref{model}), it is reasonable to consider the efficiency of our approximation space within the context of single-phase and further estimate the effect on the two-phase model. Both the two models are solved in the domain $\Omega=[0,1]\times[0,1]$.

\subsection{Single-phase flow}
In the single-phase model, we solve the pressure equation (\ref{model}) with $\lambda=1$ and postprocess the velocity field with the technique introduced in section \ref{postprocess}. For boundary condition, we set Dirichlet boundary condition $p_{D}=1$ on the left edge and $p_{D}=0$ on the right edge of the domain. Besides, we set $g_{N}=0$ on $\Gamma_{N}$, i.e. zero Neumann boundary condition for bottom and top edge. We assume there is no external force so we take $q=0$.  The size of the permeability coefficient is $100*100$ and $200*200$ for $\kappa_2$ and $\kappa_3$ accordingly. To estimate our method GMsFEM, we compare four method, the standard finite element method and GMsFEM with different combination of multiscale basis functions. For GMsFEM, we set the coarse mesh size to be $10*10$. In other words, there are $20*20$ coarse elements in the whole domain. As is shown in \ref{tab:v_error}, there is significant error decay in both cases, where online enrichment contributes much more compared to the offline enrichment. As for the notation, we use $L_z=a+b$ to denote the case where a offline basis followed by b online basis are used in each local neighborhood. In particular, with $\kappa_2$, we can see a sharp decrease from the initial case with big error to a relatively low error when we add the both offline and online basis to the case $2+1$, which is even slightly lower than case $5+0$. In other words, the information contained in online basis functions results in bigger help than that in offline basis. This is due to the global construction of online basis functions while the offline space is constructed locally. As to the other case with $\kappa_3$, we can see the decay is less pronounced, however, there is still evident improvement in the accuracy. It is similar here that we can use online basis functions to obtain satisfying results with smaller dimension of multiscale space.

To better present the approximation of the velocity field, which is actually used in the further two-phase flow, we plot figure \ref{v2} and \ref{v3} for $\kappa_2$ and $\kappa_3$. In both cases, we show the horizontal and vertical components of velocity. Since the velocity is highly related with the permeability, we exhibit the velocity field under some region with big contrast in permeabilty coefficients, which are also shown as background. We can observe significant dismatch between the initial case and reference while in the case $L_z=2+1$ and $5+0$, the accuracy improvement is apparent especially in the selected region, where the permeability changes rapidly. Specifically, in \ref{v2}, there is a few flows with opposite direction for $L_z=1+0$ compared with the reference, while the difference is less noticeable in the latter two cases.
\begin{figure}[!htbp]
    \includegraphics[width=5in]{graphs/v_k2.png}\caption{$\kappa_2$. Velocity computed using four methods. The first two columns (horizontal the first and vertical the second ) exhibits velocity profile on the whole domain with the reference on the first row and last three corresponding to $L_z=1$, $L_z=5$ and $L_z=2+1$ respectively. The last column shows the velocity  with permeability in the selected region.}
    \label{v2}
\end{figure}
\begin{figure}[!htbp]
    \includegraphics[width=5in]{graphs/v_k3.png}\caption{$\kappa_3$. Velocity computed using four methods. The first two columns (horizontal the first and vertical the second ) exhibits velocity profile on the whole domain with the reference on the first row and last three corresponding to $L_z=1$, $L_z=5$ and $L_z=2+1$ respectively. The last column shows the velocity  with permeability in the selected region.}
    \label{v3}
\end{figure}

\begin{table}[]
    \centering
    \begin{tabular}{c|c|c}
    \hline\hline
                  & $\kappa_2$ & $\kappa_3$ \\
                  \hline
       $L_z$=1+0  &  1.30  & 0.21\\
       $L_z$=5+0  &  0.16  & 0.09\\
       $L_z$=2+1  &  0.14  & 0.08\\
    \end{tabular}
    \caption{Relative error for the velocity field. The $L^2$ error of the velocity is computed for different choices of GMsFEM compared with reference solution obtained in standard finite element method.}
    \label{tab:v_error}
\end{table}
\subsection{Two-phase flow}
In solving \ref{transport}, we use the quadratic relative permeability curves $\kappa_{rw}=S^2$ and $\kappa_{ro}=(1-S)^2$, along with $\mu_w=1$ and $\mu_o=5$ for the water viscosities. The domain $\Omega=[0,1]\times[0,1]$. For the initial condition, the value at the left edge is set as $S = 1$ and we assume $S(x, 0) = 0$ elsewhere.  In practical, we construct the multiscale basis functions within the context of single-phase flow model and apply the resulted approximation space to the interested two-phase problem without updating basis function. In other words, we can precompute the bases as preparation before the simulation, which is efficient compared to the case when we need to repeat the computation for different cases.

For better visual comparison, we present the saturation for three different time levels in figure \ref{k1},\ref{k2},\ref{k3} and \ref{k4}. From figure \ref{k2} and figure \ref{k3}, significant difference from reference saturation can observed when $L_z=1$ while $L_z=2+1$ are relatively indistinguishable from reference. For case $\kappa_1$ and $\kappa_4$, the improvements are less significant yet pronounced since there are few noticeable differences between last row and reference row.
In other words, online basis functions efficiently improves accuracy compared with offline case. In figure \ref{k4}, $L_z=2+1$ is a better approximation of reference even than $L_z=8+0$. As we can verify it from figure \ref{error}, there are sharp decreases by enriching multiscale space from intial state especially with $\kappa_2$ and $\kappa_3$ compared to the other two cases, which is consistent with the previous dynamics of saturation.
Furthermore, the relative errors are improved by increasing the number of $L_z$ up to a certain threshold, and then the reduction is minimal as more functions are added. In particular, in figure (d), when we double the number of offline basis functions from the very beginning, i.e. single in each local neighborhood, the error reduction is evident however improvement is indistinguishable from $2+0$ to $8+0$, which means very limited reduction can be expected by further increasing offline basis functions. At the same time, adding a few online basis functions will notably increase the accuracy since the error in $L_z=2+2$ is even lower than the case $L_z=8$, which shows the power of incorporating global information inherited in the online basis functions. For $\kappa_1$ and $\kappa_3$, the reduction resulted by enrichment is relatively steady while in $\kappa_2$, it is easier to reach a threshold. This is due to higher heterogeneity in $\kappa_1$ compared with $\kappa_2$. However, all the above four cases combined with the single-phase case share the same conclusion that online enrichment offer us better accuracy with relatively lower cost. Therefore, it is efficient to compute residual-driven online basis functions for the sake of increasing accuracy. \\
\begin{figure}[!htbp]
    \includegraphics[width=5in]{graphs/s/s_k1.png}

    \caption{$\kappa_1$. The reference saturation is shown on the top row at three different time levels. The second and third rows are $L_2=1$, $L_z=5$ respectively.The last row is using 2 offline basis and enriched by one online basis for each coarse neighborhood, which is denoted by $L_z=2+1$.}
    \label{k1}
\end{figure}
\begin{figure}[!htbp]
    \includegraphics[width=5in]{graphs/s/s_k2.png}\caption{$\kappa_2$. The reference saturation is shown on the top row at three different time levels. The second through fourth rows are $L_z=1$, $L_z=5$ and $L_z=2+1$ respectively.}
    \label{k2}
\end{figure}
\begin{figure}[!htbp]
    \includegraphics[width=5in]{graphs/s/s_k3.png}\caption{$\kappa_3$. The reference saturation is shown on the top row at three different time levels. The second through fourth rows are $L_z=1$, $L_z=5$ and $L_z=2+1$ respectively.}
    \label{k3}
\end{figure}
\begin{figure}[!htbp]
    \includegraphics[width=5in]{graphs/s/s_k4.png}\caption{$\kappa_4$. The reference saturation is shown on the top row at three different time levels. The second through fourth rows are $L_z=3$, $L_z=8$ and $L_z=2+1$ respectively.}
    \label{k4}
\end{figure}
\begin{figure}[!htbp]
\subfigure[$\kappa_1$]
{\includegraphics[width=0.48\textwidth]{graphs/error/error_k1.png}}
\subfigure[$\kappa_2$]
{\includegraphics[width=0.48\textwidth]{graphs/error/error_k2.png}}
\subfigure[$\kappa_3$]
{\includegraphics[width=0.48\textwidth]{graphs/error/error_k3.png}}
\subfigure[$\kappa_4$]
{\includegraphics[width=0.48\textwidth]{graphs/error/error_k4.png}}
\caption{Comparison of the $L^2$ error of the saturation for $\kappa_1$ (upleft) ,$\kappa_2$ (upright) ,$\kappa_3$ (downleft) ,$\kappa_4$ (downright)as a function of time .}\label{error}
\end{figure}
\section{Conclusion}
In this paper, we consider a conservation GMsFEM for treating the coupled pressure-convection-diffusion system in the context of the two-phase flow model. An advantage of the proposed method is that local conservative and accurate velocity field can be obtained, which means we combined two main procedures, postprocessing and online enrichment. The effect can be verified in the numerical results. In the future, we can work on more computational-efficient ways to achieve the goal.

\bibliographystyle{IEEEtran}
\bibliography{reference}
\end{document}