Is there any relation between weak and strong forces, similar to electric and magnetic forces? Is it possible to unify the strong, weak, electric and magnetic field just by Maxwellian-type equations?
(Maxwell by adding a small change - unified electric and magnetic field, then Einstein's equations - use it to create special and general theory of relativity, now maybe all we need is a little more change to unify all fields,)
When $B$ and $E$ are known magnetic and electric fields, then $W$ and $S$ are weak and strong fields, what are their units?
 A: I'm guessing at what you're asking - ignore this answer if I've misunderstood you.
Maxwell's equations describe the classical behaviour of electromagnetism. They can only do this because the EM forces are long range so at macroscopic distances they behave classically. By contrast, the weak and strong forces are short range and cease to act over anything like the classical limit of distance. There is no classical approximation to describe the weak and strong forces, so there is no analogy to Maxwell's equation.
Above the electroweak transition the electroweak force will become long range. Whether there is some classical limit in the spirit of the Maxwell's equations is a good question and I don't know the answer. I would guess that there is no such limit for the strong force even about the EW transition because the force will still be confined.
A: The answer to the question posed in the title, itself, is no. The unifying connection, analogous to Electricity + Magnetism ⇒ Electromagnetism, that you're seeking out is not Weak Nuclear Force + Strong Nuclear Force ⇒ (???), but Weak Nuclear Force + Electromagnetism ⇒ Electroweak Force. Only the similarity of the names suggest otherwise.
Instead, the strong force ramifies from the deeper, more fundamental (and far stronger) quark or color force, the corresponding quantum field theory being Quantum Chromodynamics or QCD; while for the electroweak force, it is Quantum Flavor Dynamics or QFD.
So, there are only 2 known fundamental forces - apart from gravity - the electroweak and color forces; and room is provided by the Standard Model and the underlying quantum field theory for the possibility of a third, which requires a lengthy explanation about chirality, helicity, the triangle anomaly, and the fermion spectrum for the Standard Model in some detail, so I won't say anything more about it here, except a brief point before going on.
The answer to the question inside the question ("are there Maxwell-like equations?") is yes, but will require some setup. And I won't lay out the answer in full, since that involves delving more deeply into the Higgs.
The brief point that needs to be made is that the symmetry group SU(3) × SU(2) × U(1) involves the symmetries SU(3) for "color", SU(2) for "isospin" (which, despite its name, has no connection to spin or angular momentum other than its symmetry group is also an SU(2)) and U(1) for "hypercharge". The first two have non-Abelian gauge fields, while the last does not. The Maxwell equations for non-Abelian gauge fields are non-linear and inhomogeneous and all have source terms arising from the field itself - including magnetic source & current terms.
More importantly: the electromagnetic field is not the U(1) hypercharge field, but a superposition of the hypercharge field with one of the 3 modes of the SU(2) field. As such, its field equations are also non-linear and inhomogeneous. The most important feature of electroweak theory is that, by virtue of this non-linearity and inhomogeneity, its equations both supersede and falsify Maxwell's equations.
For the classical field theory, the situation is now that the Maxwell equations for electromagnetism are replaced by equations that are part of the gauge field equations written in Maxwell form. When we quantize electromagnetism within QFD, this is what we are now quantizing; not Maxwell's equations.
The Fields
The description to follow is generic to all gauge fields, irrespective of what its symmetry group is or what type of Lagrangian describes its dynamics (Yang-Mills or otherwise). The only assumption made is that there be an action principle for its dynamics governed by a Lagrangian.
The field and equations are basically the same as those of Maxwell, except that each component is replicated. For the symmetry group U(1) × SU(2) × SU(3), which has 12 generators, that means a 12-fold replication.
For the following, they are Lie vector-valued: the electric potential $φ$, magnetic potential $ = (A₁, A₂, A₃)$, magnetic induction $ = (B¹, B², B³)$ and electric force $ = (E₁, E₂, E₃)$. The scalar and each component of the three vectors is a Lie vector.
For the following, they are Lie co-vector valued: the electric induction $ = (D¹, D², D³)$, the magnetic force $ = (H₁, H₂, H₃)$, the charge density
$ρ$, and the current density $ = (J¹, J², J³)$, as well as the charge $e$. Likewise, that means every scalar and every component of every vector.
They arrange themselves into the following Lie vector valued forms:
$$A ≡ A₁ dx + A₂ dy + A₃ dz ‒ φ dt,$$
$$F ≡ B¹ dy∧dz + B² dz∧dx + B³ dx∧dy + (E₁ dx + E₂ dy + E₃ dz)∧dt,$$
respectively for the potential 1-form, field strength 2-form, and the following Lie co-vector valued forms:
$$G ≡ D¹ dy∧dz + D² dz∧dx + D³ dx∧dy ‒ (H₁ dy + H₂ dy + H₃ dz)∧dt,$$
$$Q ≡ ρ dx∧dy∧dz ‒ (J¹ dy∧dz + J² dz∧dx + J³ dx∧dy)∧dt.$$
respectively for the response field 2-form and source current 3-form.
In his earliest papers, before the treatise, Maxwell wrote the fields as differential forms, except he didn't combine the $(φ, ), (, ), (, ), (ρ, )$ pairs with $dt$, but kept them separate, and didn't use Grassmann algebra for the differential forms (except once in the treatise). And he did, indeed, briefly describe the possibility of gravity and electromagnetism combining with a non-trivial mixing angle into what we would now call an Abelian gauge field for U(1) × U(1). The modern-day analogue of this would arise by combining the hypercharge U(1) with the U(1) of the extra force I alluded to briefly above (the $B-L$ force, or force which would couple to baryon number minus lepton number), and would be describing the unification of the electroweak force with the $B-L$ force, rather than unification of electromagnetism with gravity.
The Algebraic Conventions
Every Lie group has a Lie algebra whose elements, "Lie vectors", may be expanded in terms of a basis $(Y₀, Y₁, Y₂, ⋯) = (Y_a)$. The fundamental operations of the Lie algebra are addition, multiplication by a scalar and the Lie bracket $[⋯,⋯]$.  The bracket operation is determined completely by its action on the basis,
$$[Y_a,Y_b] = Σ_c f_{ab}^c Y_c,$$
expressed in terms of the structure coefficients $f_{ab}^c$. In here, and the following, I will used the Einstein summation convention, which means that repeated indices are to be summed over. So the equation above can be rewritten as $[Y_a,Y_b] = f_{ab}^c Y_c$.
Consequently, the fields, their components and the forms $A$, $F$, $G$ and $Q$, and the charge $e$ all have decompositions in terms of the respective bases
$$(φ, , , ) = (φ^a Y_a, ^a Y_a, ^a Y_a, ^a Y_a),$$
$$(ρ, , , ) = (ρ_a Y^a, _a Y^a, _a Y^a, _a Y^a),$$
$$e = e_a Y^a.$$
The important point to make about this is that the two members in each of the pairs $(,)$ and $(,)$ cannot be simply equated to one another, because they're not even the same type of objects.
It is always possible to expand a Lie algebra by endowing it with an associative product operation such that $uv ‒ vu = [u,v]$; for instance by adopting a faithful matrix representation or an embedding algebra.
The dual of the Lie algebra has a basis $(Y⁰, Y¹, Y², ⋯) = (Y^a: a = 0, 1, 2, ⋯)$. The fundamental operations are contraction and a "co-adjoint bracket" with Lie vectors:
$$Y_a ˩ Y^b = δ_a^b,$$
$$[Y^c, Y_a] = f^c_{ab} Y^b$$
where $δ_a^b = 1$ if $a = b$ and $δ_a^b = 0$ if $a ≠ b$. The latter operation satisfies the identity
$$v ˩ [ω, u] = [u, v] ˩ ω.$$
For semi-simple Lie groups, the algebraic extension for their Lie algebras can be further extended to include co-vectors, such that
$$ω u ‒ u ω = [ω, u]$$
and an additional operator, a linear "trace" operator that satisfies:
$$Tr(ωu) = u ˩ ω,$$
$$Tr(ab⋯c) = Tr(b⋯ca).$$
$$Tr(u ± v) = Tr(u) ± Tr(v)$$
Thus, for instance
$$Tr(ω [u,v]) = Tr(ω (uv ‒ vu)) = Tr((ωu ‒ uω)v) = Tr([ω,u]v)$$
which is consistent with the requirement that $v ˩ [ω,u] = [u,v] ˩ ω$.
By virtue of the Lie algebra being embedded in an associative algebra, as described above, the vector operations can be nicely dovetailed into the Lie operations. Thus,
$$ ×  = ^a × ^b Y_a Y_b = ½ ^a × ^b (Y_a Y_b ‒ Y_b Y_a) = ½ f^c_{ab} ^a × ^b Y_c$$
and, similarly:
$$φ ‒ φ = f^c_{ab} φ^a ^b Y_c;$$
and for the co-vector valued forms, we can write
$$φ ‒ φ = -f^c_{ab} φ^a _c Y^b.$$
Gauge Field Equations In Maxwell Form
With these conventions, the equations for classical gauge theory, may thus be written in a form analogous to Maxwell's equations as:
1. Field-Potential Equations:
$$ = ∇ ×  +  × ,$$
$$ = ‒\frac{∂}{∂t} ‒ ∇φ + φ ‒ φ.$$
2. Bianchi Identities:
$$∇· + · ‒ · = 0,$$
$$∇× + \frac{∂}{∂t} + × + × + φ ‒ φ = .$$
3. Field-Law:
$$∇· + · ‒ · = ρ,$$
$$∇× ‒ \frac{∂}{∂t} + × + × + φ ‒ φ = .$$
4. Continuity Equation:
$$∇· + \frac{∂ρ}{∂t} + · ‒ · + ρφ ‒ φρ = · ‒ · + · ‒ ·.$$
The right hand side is always zero for field theories whose dynamics are derived from Lorentz invariant Lagrangians.
5. Force and Power Laws (first form):
$$\frac{d( + Tr(e))}{dt} = ‒∇(Tr(e(φ ‒ ·))),$$
$$\frac{d(H + Tr(eφ))}{dt} = \frac{∂(Tr(e(φ ‒ ·)))}{∂t}.$$
both partial derivatives on the right being taken, holding $e$ and $ = d/dt$ fixed.
6. Force and Power Laws (second form):
$$\frac{d}{dt} = Tr(e( + ×)),$$
$$\frac{dH}{dt} = Tr(e·),$$
$$\frac{de}{dt} = Tr((φ‒·)e ‒ e(φ‒·)),$$
the last one being called Wong's Equation, which only holds if the right-hand side of the continuity equation is 0; i.e. if $· ‒ · + · ‒ · = 0$.
Wong's equation describes the precession of the gauge charge under the influence of the gauge field. For isospin SU(2), this precession is between the isospin-up and isospin-down states, which correspond respectively to ${u,c,t,ν_e,ν_μ,ν_τ}$ and their respective partners ${d,s,b,e,μ,τ}$. (And it only applies to left helicity states for matter and right helicity states for anti-matter, which is what I didn't want to get into in detail here.)
7. Force and Power Laws (continuum form):
$$ = Tr(ρ + ×),$$
$$P = Tr(·).$$
Equations 1-4 are respectively the component forms of
$$dA + A² = F,$$
$$dF + AF - FA = 0,$$
$$dG + AG - GA = Q,$$
$$dQ + AQ + QA = FG - GF,$$
where similar efficiencies occur by combining the Grassmann algebra of the differential forms with the algebra for the Lie vectors and covectors.
The other equations (5, 6 and 7) can also be written in terms of natural operations with differential forms, but are a little more involved, and I won't get into it here. As such, they are all natural equations that do not require any background metric, and are independent of metrical and causal structure. They take on the same form, regardless of whether the underlying manifold is endowed with a 4+0 locally Euclidean metric, a 3+1 locally Minkowski metric, or even a Newton-Cartan structure for non-relativistic field theory.
The Constitutive Laws with Yang-Mills Lagrangians
8. In-vacuuo constitutive Law
$$_a = ε_{0ab} ^b,$$
$$^a = μ_0^{ab} _b,$$
where $k_{ab} = ε_{0ab} c$ and $k^{ab} = μ_0^{ab} c$ are, respectively, the components of the metric and inverse metric for the underlying Lie algebra with respect to the bases $(Y_a)$ and $(Y^a)$.
This arises, as the field equations, of the Yang-Mills action
$$S = ∫ ½ \frac{k_{ab}}{c} (^a·^b ‒ ^a·^b c²) d⁴x$$
which is the direct generalization of the Maxwell-Lorentz action
$$S = ∫ ½ ε₀ (E² ‒ B² c²) d⁴x$$
and it shows that the coefficient $k ≡ ε₀ c = \sqrt{ε₀/μ₀}$ is actually the gauge group metric for U(1), when the Maxwell-Lorentz action is written as a Yang-Mills action for U(1). It is not a coefficient that is "purely conventional" nor one that can be "ignored" or "defined away", but has physical and geometric meaning: it describes a crucial part of the very geometry of U(1), itself ― its metric!
The physical relevance of the constitutive coefficients and constitutive law is a point that Hehl has made in the context of electromagnetic theory. It contains important physics that will come back to bite you (in disguise) below, in the context of Quantum Field Theory, even if you try to remove it from the classical theory. More importantly, in non-Cartesian coordinates and in curved spacetimes, the coefficient $ε₀c$ becomes a function of the space-time metric and coordinates. As Hehl pointed out, this then describes - in a non-trivial way - constitutive property of the vacuum itself. The vacuum is no longer something that has nothing in it. The metric counts as "something" (as does the connection, but that doesn't directly come into play here), not "nothing".
It is more common to see the fields written in the theoretical literature in Gaussian form (denoted here with the subscript $()₊$):
$$(₊, φ₊, ₊, ₊) = (\sqrt{\frac{4π}{μ₀}} , \sqrt{4πε₀}φ, \sqrt{\frac{4π}{μ₀}} , \sqrt{4πε₀}),$$
$$(₊, ρ₊, ₊, ₊) = (\frac{}{\sqrt{4πε₀}}, \frac{ρ}{\sqrt{4πε₀}}, \sqrt{\frac{4π}{ε₀}} , \sqrt{4πμ₀}),$$
$$(e₊) = (\frac{e}{\sqrt{4πε₀}}),$$
$$(μ₊, ε₊) = (\frac{μ}{μ₀}, \frac{ε}{ε₀}).$$
with or without the extra $4π$'s.
This has the unfortunate effect of muddying the picture, by mixing in components of the gauge group metric with the fields themselves. Scale-dependency of the metric under renormalization is then mixed up with and confused with scale dependency of the field components and charge - one of the biggest errors and misconceptions propagated in the literature (although, there are a few places in the literature where it's been called out with both regularization and renormalization described in terms of an effect on the constitutive coefficients). It's actually the metric which is scale-dependent, rather than the field and charge.
The way you can tell that the "scale dependency" of renormalization theory actually pertains to the gauge group metric, itself, is that the renormalization for the potential $A$ and field strength $F$ goes like $A₊ ⇒ √Z₃ A₊$, $F₊ ⇒ √Z₃ F₊$ when the fields are written under Gaussian convention, and if you closely examine both this and the Lagrangian, you will see that the Lagrangian goes like $½ ε₀ (E² ‒ B²c²) ⇒ ½ Z₃ ε₀ (E² ‒ B²c²)$, once the Gaussian scaling is removed. Likewise, the Gaussian scaling can all be written with $√ε₀$ in the numerator, as
$$(₊, φ₊, ₊, ₊) = \sqrt{4πε₀} (c, φ, c, ) ⇒ \sqrt{4πε₀Z₃} (c, φ, c, ) = \sqrt{Z₃} (₊, φ₊, ₊, ₊).$$
So, it's actually $ε₀$ that $Z₃$ is going with, rather than with $A$, $F$ or its components!
Similarly, for gauge theory, the coefficient $Z₃$ belongs with the gauge group metric $k_{ab}$, which is what $ε₀c$ corresponds to. To rescale the fields in a way analogous to the Gaussian convention would require first taking the square root of the metric (i.e. a singular value decomposition of the metric). So, it is also normally assumed that the gauge group metric is positive or negative definite.
The Constitutive Laws in General with Lorentz-Invariant Lagrangians
The constitutive laws posed under item 8 are those for the Yang-Mills action. More generally, for a gauge theory that arises from an action given by
$$S = ∫ (,φ,,) d⁴x$$
with Lagrangian density $$, the constitutive laws may be written:
$$_a = \frac{∂}{∂^a}, _a = -\frac{∂}{∂^a}, _a = \frac{∂}{∂^a}, ρ_a = -\frac{∂}{∂φ^a}$$
Applying the variational
$$δ = Tr(δ· ‒ δ· + δ· ‒ δφ ρ)$$
and using the Field-Potential Equations to reduce $δ$, $δ$ to $δ$, $δφ$
$$δ = Tr((-∇δφ ‒ \frac{∂(δ)}{∂t} + δ(φ ‒ φ))· ‒ (∇×(δ) + δ(×))· + δ· ‒ δφ ρ)$$
integrate by parts to get:
$$δ = Tr(∇(-δφ· ‒ δ×) ‒ \frac{∂(δ·)}{∂t} + δφ(∇· + · ‒ · ‒ ρ) ‒ δ·(∇× ‒ \frac{∂}{∂t} + × + × + φ ‒ φ ‒ ))$$
and obtain, as a result, the field equations.
If the Lagrangian is Lorentz-dependent, then the field-dependent part of the Lagrangian would be a function only of its Lorentz invariants
$$ℑ^{ab} ≡ ½ (^a·^b ‒ ^a·^b c²),$$
$$^{ab} ≡ ½ (^a·^b + ^a·^b),$$
with the resulting constitutive law:
$$_a = ε_{ab} ^b + θ_{ab} ^b,$$
$$_a = ε_{ab} c² ^b ‒ θ_{ab} ^b,$$
which includes an axial version $θ_{ab}$ of permittivity, where the coefficients are functions of $ℑ$ and $$
$$ε_{ab}(ℑ,) = \frac{∂}{∂ℑ^{ab}}, θ_{ab}(ℑ,) = \frac{∂}{∂^{ab}}$$
such that
$$\frac{∂ε_{ab}}{∂ℑ^{cd}} = \frac{∂ε_{cd}}{∂ℑ^{ab}},
\frac{∂ε_{ab}}{∂^{cd}} = \frac{∂θ_{cd}}{∂ℑ^{ab}},
\frac{∂θ_{ab}}{∂^{cd}} = \frac{∂θ_{cd}}{∂^{ab}}.$$
They need not be constant. In a quantized theory, they would be modelled as "constant" $ε₀$, $θ₀$ that take on different appearances at different scales ― i.e. the very essence of "renormalization", itself ― and would produce a corresponding renormalization group flow for the fields and charge, when they are written in Gaussian form.
A consequence of these constitutive laws is that the right-hand side of the continuity equation is 0:
$$· ‒ · + · ‒ · = 0.$$
Therefore, the continuity equation reduces to
$$∇· + \frac{∂ρ}{∂t} + · ‒ · + ρφ ‒ φρ = 0,$$
and Wong's equation (in 6) holds.
The $θ$ coefficients are defined only up to an additive constant, since the field equations are symmetric under the transform
$$(_a, _a) → (_a ‒ θ_{0ab} ^b, _a + θ_{0ab} ^b),$$
which is a residue and reduced version of complexion symmetry.
For electromagnetism, the conditions $(ℑ,) = (0,0)$ define the null field - i.e. the field associated with pure radiation. So, even for general Lagrangians, we can define the null field values of the constitutive coefficients (suppressing Lie algebra indices):
$$ε₀ ≡ ε(0,0), k₀ = ε₀ c, θ₀ ≡ θ(0,0) = 0.$$
the latter can be set to 0, by using the above transform. Correspondingly, by Taylor's Theorem (that is: the exact form of Taylor's Theorem - with remainder), it follows that all Lagrangian densities that involve gauge fields - if the Lagrangian desnity possesses Lorentz symmetry - can be written as
$$(ℑ,) = (0,0) + k₀ ℑ + ½ (β₀ ℑ ℑ + β₁ (ℑ  +  ℑ) + β₂  ).$$
For weak fields, this reduces to
$$(ℑ,) ≃ _M + _{YM}$$
where
$$_M ≡ (0,0), _{YM} = k₀ ℑ$$
is the Lagrangian density for matter and for all the other fields and $_{YM}$ is the Yang-Mills Lagrangian itself.
The Constitutive Laws for Isotropic Media
If the Lagrangian is only isotropic (as would be suitable for isotropic media), rather than also boost-invariant, then the field-dependent part would be a function of its scalar invariants. In that case, it would actually be better to revert to a treatment closer to Maxwell, in spirit, by using the Routhian density $ℜ =  + Tr(·)$, in place of $$, with the variational
$$δℜ = Tr(δ· + δ· + ⋯)$$
that treats the "induction" fields $$, $$ as derived and the "force" fields $$, $$ as fundamental. The corresponding scalar invariants are
$$ℑ^{ab} = ^a·^b, ^a_b = ½ ^a·_b, _{ab} = _a·_b$$
with coefficients given by
$$κ_{ab} = \frac{∂ℜ}{∂ℑ^{ab}}, λ_a^b = \frac{∂ℜ}{∂^a_b}, μ^{ab} = \frac{∂ℜ}{∂_{ab}}$$
and constitutive law
$$_a = κ_{ab} ^b + λ_a^b _b,$$
$$^a = λ_b^a ^b + μ^{ab} _b.$$
The symbols for the dielectric coefficient $κ$ and permeability $μ$ actually go all the way back to Maxwell, who devised the names by analogy to the spring coefficient $k$ and mass $m$. The axial coefficient $λ$ was not part of anyone's exposition of classical theory, but is included in the list of possibilities as is $θ$. Both are pseudo-scalars, and if they are 0, then both $κ$ and $ε$ coincide: $θ_{ab} = 0 ⇔ λ_a^b = 0 ⇒ κ_{ab} = ε_{ab}$.
When in vacuo, they are related to the first set of constitutive coefficients by:
$$ε_{ab} = κ_{ab} ‒ λ_a^c μ⁻¹_{cd} λ_b^d,$$
$$θ_{ab} = λ_a^c μ⁻¹_{cb},$$
$$ε_{ab} c² = μ⁻¹_{ab}.$$
This relation is only consistent if $μ$ is non-singular; i.e. if the Hessian $\frac{∂²ℜ}{∂_a ∂_b}$ is non-singular. But that's already the precondition for being able to convert between the Routhian density $ℜ$ and Lagrangian density $$.
Isotropy is a frame-dependent condition, in general. Therefore, these relations can only be posed in one frame. In Maxwell's time and the early 20th century, this was called the stationary frame. For other frames, they would take the form of
The Generalized Maxwell-Minkowski Relations:
$$ ‒ α  ×  = λ ( +  × ) + μ ( ‒  × ),$$
$$ + α  ×  = κ ( +  × ) + λ ( ‒  × ).$$
For general values of $α$, these are the equations suitable for a geometry which has the following invariants:
$$dt² ‒ α(dx² + dy² + dz²)$$
$$((∂/∂x)² + (∂/∂y)² + (∂/∂z)²) ‒ α (∂/∂t)²$$
It specializes to the Maxwell-Minkowski relations when $α > 0$, with light speed given by $c = \sqrt{1/α}$; to their non-relativistic version when $α = 0$; and to 4D Euclidean form $α < 0$, which could be used for Euclidean field theory.
In the Maxwell-Minkowski relations, invariance under boosts are (generally) broken, and the velocity $$ marks out and identifies the distinguished frame in which the constitutive relations reduce to isotropic form $ = λ + μ$, $ = κ + λ$.
All of this generalizes accordingly to gauge fields, with the coefficients $κ$, $λ$ and $μ$ replaced, respectively by $κ_{ab}$, $λ^a_b$ and $μ^{ab}$:
$$^a ‒ α  × ^a = λ^a_b (^b +  × ^b) + μ^{ab} (_b ‒  × _b),$$
$$_a + α  × _a = κ_{ab} (^b +  × ^b) + λ^b_a (_b ‒  × _b).$$
The product matrix $εμ = κμ ‒ λλ^T$ would then have a singular value decomposition that, generalizing the relation $εμ = 1/V²$, which yields wave speed $V$ in the medium, would produce eigenvalues that give different wave speeds for different eigenmodes.
I've never seen anyone do anything with classical or quantum gauge fields for moving media, so this may all be new. But it is the obvious generalization of the classical theory of moving media for electromagnetic theory and it may be possible to tie it in directly with the Effective Lagrangian approach used in Quantum Field Theory and actually give some precise meaning to and quantify the notions of "screening" or "anti-screening" vacuua and other forms of effective media.
