Increasing a potential causes increase in energy levels Suppose a potential $V(x)$, and suppose a bound particle so the allowed energy levels are discrete. Suppose a second potential $\widetilde{V}(x)$ such that $\widetilde{V}(x) \geq V(x)$ for all $x$ (suppose the potentials are relevant only for some interval). Does this necessarily imply that the eigenvalues of the Hamiltonian with $\widetilde{V}(x)$ will be at least those of the Hamiltonian with $V(x)$? That is, if $E_n$ are the energy levels of the first Hamiltonian and $\widetilde{E}_n$ are the energy levels of the second Hamiltonian, is $E_n\leq \widetilde{E}_n$ for all $n$? How can one show this?
 A: You can show this by using perturbation theory (only for suitable small changes in the potential).
When you assume, that $\tilde{V}(x) = V(x) + c$ with $c > 0$, then you can write your problem als perturbation:
If the unperturbated hamiltonian $\hat{\mathrm H}$ has eigenstates $ | \Psi_n \rangle $ with discrete energies, then perturbation-theory states that changing the hamiltonian by a little term $\hat{\mathrm V}_\textrm{perturbation}$ will change the eigenvalues $E_n$ by:
$$
\Delta E_n = \left\langle \Psi_n \left| \hat{\mathrm V}_\textrm{perturbation}\right| \Psi_n \right\rangle
$$
This is valid if you neglect terms of higher Order. 
There is one think to watch out regarding perturbation theory: If your energy eigenvalues are degenerated, then the perturbation term has to be diagonal in the subspace that is spanned by the degenerated states. 
In our case, $ \hat{\mathrm V}_\textrm{perturbation}= c$ is just a multiplication, so:
$$
\langle \Psi_m | c | \Psi_n\rangle = c \langle \Psi_m | \Psi_n\rangle = c ~\delta_{nm}
$$
$\hat{\mathrm V}_\textrm{perturbation}$ is diagonal in any subspace, and we can make use of perturbation theory. You then calculate the energy-shift just by
$$
\Delta E_n = \langle \Psi_n | c| \Psi_n \rangle = c  \langle \Psi_n | \Psi_n \rangle = c >0
$$
So if you increase the potential by a constant, then the energy eigenvalues will just shift by that constant. 
Edit: One can expand the proof for perturbations that vary with time: Let the change in the potential be $\delta V(x)$ (which now depends on $x$), then you can still calculate the energy-shift by using perturbation theory. In whatever subspace that is formed by degenerate states, you can find a Basis $|\tilde{\Psi}_n \rangle$ for which $\delta V(x)$ is an orthogonal Operator. 
In this Basis you then calculate the energy-shift like described above:
$$
\Delta E_n = \langle \tilde{\Psi}_n | \delta V( \hat{x})| \tilde{\Psi}_n = \int dx |\tilde{\Psi}(x)|^2 \delta V(x)   > 0
$$
Since $\delta V(x) > 0$. Those are now energyshifts for eigenstates of your "old" hamiltonian. However, those eigenstates are not necessarrily the eigenstates that you started with.  
A: *

*Let there be given a selfadjoint$^1$ operator $H^{(0)}$ and a (semi)positive operator $V\geq 0$ on a Hilbert space ${\cal H}$. Let the basis of normalized eigenvectors for $H^{(0)}$ be $(|i^{(0)}\rangle)_{i\in I}$ with corresponding eigenvalues $(E^{(0)}_i)_{i\in I}$ ordered such that 
$$ \forall i,j ~\in~I:\quad i~\leq~j\quad\Rightarrow \quad E^{(0)}_i~\leq~E^{(0)}_j.\tag{1}$$ Similarly, let the basis of normalized eigenvectors for $H^{(1)}:=H^{(0)}+V$ be $(|i^{(1)}\rangle)_{i\in I}$ with corresponding ordered eigenvalues $(E^{(1)}_i)_{i\in I}$.

*Non-degenerate$^2$ perturbation theory yields the following construction: Define a one-parameter family of selfadjoint operators $$H(t)~:=~ H^{(0)}+tV, \qquad t~\in~[0,1].\tag{2}$$ Consider the following initial value problem of coupled 1st order differential equations$^2$
$$ \frac{d|i(t)\rangle}{dt}~~=~~\sum_{j\in I\backslash \{i\}} \frac{\langle j(t)| V | i(t) \rangle}{ E_i(t)-E_j(t)}|j(t)\rangle ~~=~~ \sum_{j\in I}A_{ij}(t)|j(t)\rangle , \tag{3}$$
$$ A_{ij}(t)~~:=~~\left\{\begin{array}{ccl} \frac{\langle j(t)| V | i(t) \rangle}{ E_i(t)-E_j(t)} &\text{if}& i~\neq~ j \cr
0 &\text{if}& i~=~ j \end{array} \right\}~~=~~-A^{\ast}_{ji}(t),  \tag{4}$$
$$ \frac{dE_i(t)}{dt}~~=~~\langle i(t)| V | i(t) \rangle~~\geq~~0, \tag{5}$$
$$|i(t\!=\!0)\rangle~~=~~|i^{(0)}\rangle,\qquad E_i(t\!=\!0)~~=~~E^{(0)}_i,\qquad i\in I. \tag{6}$$
We conclude from eq. (5) that the spectrum of $H(t\!=\!1)\equiv H^{(1)}$ is increased compared to the spectrum of $H(t\!=\!0)\equiv H^{(0)}$, i.e. 
$$ \forall i~\in~I:\quad E^{(0)}_i ~\leq~ E^{(1)}_i,\tag{7}$$ 
as OP wanted to know. Note that the anti-Hermitian property (4) implies that the basis $(|i(t)\rangle)_{i\in I}$ is normalized. 

*In case of degeneracies and level-crossings, the basis $(|i(t)\rangle)_{i\in I}$ is no longer a well-defined/continuous function of $t$. Nevertheless, one may argue that the ordered eigenvalues $E_i(t)$ are still non-decreasing as a function of $t$, i.e. weakly increasing.

*Alternatively, the increase in the ground state energy can be independently & non-perturbatively deduced from the variational method:
$$ E_0^{(0)} ~~\leq~~ \langle 0^{(1)}| H^{(0)} | 0^{(1)} \rangle
~~=~~ \langle 0^{(1)}| \left(H^{(1)}\!-\! V \right)| 0^{(1)} \rangle~~\leq~~ E^{(1)}_0. \tag{8}$$

*A generalization of the above variational method (8), based on the Schur-Horn theorem, leads to the following tower of weaker inequalities:
$$ \forall j ~\in~I:\quad \sum_{i=0}^j E^{(0)}_i ~\leq~\sum_{i=0}^j E^{(1)}_i.\tag{9}$$ 
(To prove ineq. (9) work in a basis where $H^{(1)}$ is diagonal, and then apply the Schur-Horn inequality.)

$^1$We will ignore subtleties with unbounded operators, domains, selfadjoint extensions, etc., in this answer.
$^2$ We assume in Section 2 that the spectrum for $H^{(0)}$ is non-degenerate and that no level-crossings occur.
A: Consider $\tilde{V}(x)=V(x)+\Delta(x)$ where $\Delta(x)>0,\forall x$. Within 1st order perturbation, 
$E_n=E_n^0+\langle \psi_n|\Delta|\psi_n\rangle=E_n^0+\iint dx_1dx_2\psi_n^*(x_1)\delta(x_1-x_2)\Delta(x_2)\psi_n(x_2)=E_n^0+\int dx|\psi_n(x)|^2\Delta(x)>E_n^0$
A: Generically, the answer is "no".  The Zeeman effect is the splitting of degenerate spectral lines in the presence of a static magnetic field.  As the field strength increases, some lines move to higher energies and some lines move to lower energies.  Example of the splitting of the $5s$ orbitals of Rubidium:

(Graph created by: Danski14.  Image used under Creative Commons Attribution-Share Alike 3.0 Unported license)
Note that the levels don't even move monotonically:  The (2,-1) and (1,-1) levels move towards zero frequency shift for small field and then away for large fields.
Similarly, the Stark effect (shifting and splitting of spectral lines by an imposed electrical field) can shift to higher and lower energies.  (Structurally similar shifts can be seen on the linked page as those shown above for the Zeeman effect.)
EDIT:
At least one commenter seems confused about what these words mean.  The Zeeman splitting occurs under a uniform field.  A uniform field has no potential, so makes no contribution to $\tilde{V}(x)$.  In detail, using the OP's potential, $$ V(x) + \text{[uniform magnetic field]} = V(x) + 0 = \tilde{V}(x)  \text{.}  $$  Consequently, Zeeman splitting occurs without changing the potential at all.  (This is admitted by the OP's use of $V(x) \leq \tilde{V}(x)$.)
A magnetic field is a vector field, so does not have spin dependence.  That is, the uniform magnetic field is not written as a function of position and spin; it is only a function of position.  The Hamiltonian contains all of the interaction terms.  If there is coupling between electrons and the magnetic field, this coupling appears in the Hamiltonian.  And in a Hamiltonian that has such a term, the spin of an electron does contribute to the sign of the effect of the electron-magnetic interaction.  The OP gives no details about his/her Hamiltonian.  If a commenter has questions about spin-dependence in the OP's Hamiltonian, the only profitable choice is to inquire with the OP.
It is perhaps unfamiliar to imagine a uniform electric field (rather than an electric field of uniform gradient).  Such a field is approximately produced in doped semiconductors.  Far from a p-n junction, the Stark effect deforms the orbital/band structure of the materials, differently depending on the dopant.  Since we are far from the junction, the net electric field is induced by the dopant (and depletion physics) and is as uniform as your model allows.  Near the junction, the local field value is intermediate and the Stark deformation of the energy levels/bands is intermediate (and a bit nonlinear since now the gradient of the electric field is not necessarily small).
So, without changing the numerical value of $V(x)$ at all, binding energies can can be shifted positively or negatively.  This is observed in real atoms.
