Let me discuss the mathematical precise derivation of the Yang-Mills equation in full generality. For this, let us first of all fix some notation:
- Let $P$ be a principal bundle over some (compact, oriented) pseudo-Riemannian manifold $(\mathcal{M},g)$ (spacetime) with structure group given by a compact (and finite-dimensional) Lie group $G$. The corresponding Lie algebra is denoted by $\mathfrak{g}$.
- Furthermore, let $\langle\cdot,\cdot\rangle_{\mathfrak{g}}$ be an $\mathrm{Ad}$-invariant inner product on $\mathfrak{g}$, or more generally, a non-degenerate $\mathrm{Ad}$-invariant and symmetric bilinaer form. For example, if $G$ is simple, then this is usually nothing else then a (negative multiple) of the Killing form.
Now, in order to formulate the action of Yang-Mills theory, we first of all take a connection $1$-form $A\in\Omega^{1}(P,\mathfrak{g})$, which corresponds to a gauge field in physics terminology. The corresponding curvature $F^{A}\in\Omega^{2}(P,\mathfrak{g})$ is defined by
$$F^{A}=\mathrm{d}A+\frac{1}{2}[A\wedge A].$$
Since we want to define an integral over the spacetime $\mathcal{M}$ and not $P$, we need to translate the curvature to a field defined on $\mathcal{M}$. This can be done in the following way: It is a general mathematical fact that there is the following isomorphism:
$$\Omega^{k}_{\mathrm{hor}}(P,\mathfrak{g})^{\mathrm{Ad}}\cong\Omega^{k}(\mathcal{M},\mathrm{Ad}(P)),$$
where $\Omega^{k}_{\mathrm{hor}}(P,\mathfrak{g})^{\mathrm{Ad}}$ is some subset of $\Omega^{k}(P,\mathfrak{g})$ of forms satisfying extra properties (which are fulfilled by $F^{A}$) and where $\mathrm{Ad}(P)=P\times_{\mathrm{Ad}}\mathfrak{g}$ denotes the "adjoint bundle", which is a vector bundle defined on some certain quotient of $P\times\mathfrak{g}$. (In mathematics, this is a particular case of so-called "associated vector bundles", which can be defined for every principal bundle and every representation $(\rho,V)$ on $G$). As a second ingredient, we have to define an inner product
$$\langle\cdot,\cdot\rangle_{\mathrm{Ad}(P)}:\Omega^{k}(\mathcal{M},\mathrm{Ad}(P))\times\Omega^{k}(\mathcal{M},\mathrm{Ad}(P))\to C^{\infty}(\mathcal{M}).$$
This can be done in the obvious way: Take $\alpha,\beta\in\Omega^{k}(\mathcal{M},\mathrm{Ad}(P))$. If we take some local frame $\{e_{a}\}_{a=1}^{\mathrm{dim}(G)}\subset\Gamma(U,\mathrm{Ad}(P))$, which is a family of local sections of $\mathrm{Ad}(P)$, such that they form at each point a basis, then we can write
$$\alpha\vert_{U}=\sum_{a=1}^{\mathrm{dim}(G)}\alpha^{a}\otimes e_{a}\hspace{1cm}\text{and}\hspace{1cm}\beta\vert_{U}=\sum_{a=1}^{\mathrm{dim}(G)}\beta^{a}\otimes e_{a},$$
where $\alpha^{a},\beta^{a}\in\Omega^{k}(U)$ are real-valued forms. Using this, we define the iner product by
$$\langle\alpha,\beta\rangle_{\mathrm{Ad}(P)}\vert_{U}:=\sum_{a,b=1}^{\mathrm{dim}(G)}\langle\alpha^{a},\beta^{b}\rangle\langle e_{a},e_{b}\rangle_{\mathfrak{g}},$$
where $\langle\alpha^{a},\beta^{b}\rangle$ denotes just the usual inner product of real-valued forms and where $\langle e_{a},e_{b}\rangle_{\mathfrak{g}}$ has to be understood point-wise. After all this preliminaries, we define the Yang-Mills action via
$$\mathcal{S}_{\mathrm{YM}}[A]:=\int_{\mathcal{M}}\Vert F^{A}_{\mathcal{M}}\Vert^{2}_{\mathrm{Ad}(P)}\,\mathrm{d}\mathrm{vol}_{g}$$
where $F_{\mathcal{M}}^{A}$ is the curvature viewed as an element of $\Omega^{2}(\mathcal{M},\mathrm{Ad}(P))$ via the above isomorphism and where $\mathrm{d}\mathrm{vol}_{g}$ denotes the usual measure on pseudo-Riemannian manifolds. In order to derive the equations of motion, we first of all observe that
$$F^{A+t\alpha}=F^{A}+t(\mathrm{d}\alpha+[A\wedge\alpha])+\frac{t^{2}}{2}[\alpha\wedge\alpha]$$
for all $\alpha\in\Omega^{1}_{\mathrm{hor}}(P,\mathfrak{g})^{\mathrm{Ad}(P)}$. As a consequence, we have that
$$F^{A+t\alpha}_{\mathcal{M}}=F_{\mathcal{M}}^{A}+t(\mathrm{d}_{A}\alpha_{\mathcal{M}})+\mathcal{O}(t^{2}),$$ where $\alpha_{\mathcal{M}}\in\Omega^{1}(\mathcal{M},\mathrm{Ad}(P))$ corresponds to the form $\alpha$ via the isomorphism explained above. As a last ingredient, we need the following fact:
$$\int_{\mathcal{M}}\langle\mathrm{d}_{A}\alpha,\beta\rangle_{\mathrm{Ad}(P)}\,\mathrm{vol}_{g}=\int_{\mathcal{M}}\langle\alpha,\mathrm{d}_{A}^{\ast}\beta\rangle_{\mathrm{Ad}(P)}\,\mathrm{vol}_{g}$$
for all $\alpha\in\Omega^{k}(\mathcal{M},\mathrm{Ad}(P))$ and for all $\beta\in\Omega^{k+1}(\mathcal{M},\mathrm{Ad}(P))$, where $$\mathrm{d}^{\ast}_{A}:\Omega^{k+1}(\mathcal{M},\mathrm{Ad}(P))\to\Omega^{k}(\mathcal{M},\mathrm{Ad}(P))$$ denotes the codifferential. This formula is basically an extension of the well-knwon theorem for real-valued forms, that states that the exterior derivative is formally self-adjoint to the codifferential with respect to a suitable $L^{2}$-inner product, to bundle-valued differential forms. Using this, we get that
$$\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{S}_{\mathrm{YM}}[A+t\alpha]=2\int_{\mathcal{M}}\,\langle\mathrm{d}_{A}\alpha_{M},F_{\mathcal{M}}^{A}\rangle_{\mathrm{Ad}(P)}\,\mathrm{vol}_{g}=2\int_{\mathcal{M}}\,\langle\alpha_{M},\mathrm{d}_{A}^{\ast}F_{\mathcal{M}}^{A}\rangle_{\mathrm{Ad}(P)}\,\mathrm{vol}_{g}\stackrel{!}{=}0.$$
Hence, by non-degeneracy of the inner product, which extends to the inner product $\langle\cdot,\cdot\rangle_{\mathrm{Ad}(P)}$, we get that
$$\mathrm{d}^{\ast}_{A}F_{\mathcal{M}}^{A}=0.$$
Using the definition of the codifferential, this is equivalent to say that
$$\mathrm{d}_{A}\ast F_{\mathcal{M}}^{A}=0,$$
which are the well-known Yang-Mills equation.
As a short comment: It is a general fact that for every connection $1$-form $A$ it holds that
$$\mathrm{d}_{A}F_{\mathcal{M}}^{A}=0.$$
This is a particular form of the "Bianchi identity". This equation together with the Yang-Mills equation can be viewed as a generalization of the well-known Maxwell's equations for electrodynamics.
Best Answer
Wu and Yang (1968) found a static solution to the sourceless SU(2) Yang-Mills equations, (please, see the following two relatively recent articles containing a rather detailed description of the solution: Marinho, Oliveira, Carlson, Frederico and Ngome The solution constitutes of a generalization of the Abelian Dirac monopole. The vector potential is given by:
$A_i^a=g \epsilon_{iaj}x^j\frac{f(r)}{r^2}$
where $f(r)$ satisfies a nonlinear radial equation (The Wu-Yang equation) obtained from the substitution this ansatz into the Yang-Mills field equations.The Wu-Yang monopole has a singularity at the origin, in which the magnetic energy density diverges. The first article contains references to phenomenological works involving the Wu-Yang monopole.