Felix Winterhalter
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
        variances = SymbolGen("variance") if params['fluctuating'] is True else params['fluctuating']
        correction = fluctuation_correction(lb_method, random_symbol(collision_rule.subexpressions, dim=lb_method.dim),
                                            variances)
+
        for i, corr in enumerate(correction):
            collision_rule.main_assignments[i] = Assignment(collision_rule.main_assignments[i].lhs,
                                                            collision_rule.main_assignments[i].rhs + corr)

+        variances = SymbolGen("variance") if params['fluctuating'] is True else params['fluctuating']
+        temperature = sp.Symbol("fluctuating_temperature") if params["temperature"] is None else params["temperature"]
+        variance_equations = fluctuating_variance_equations(method=lb_method,
+                                                            temperature=temperature,
+                                                            c_s_sq=params["c_s_sq"],
+                                                            variances=variances)
--- a/lbmpy/creationfunctions.py
+++ b/lbmpy/creationfunctions.py
    collision_rule = simplification(collision_rule)

    if params['fluctuating']:
-        variances = SymbolGen("variance") if params['fluctuating'] is True else params['fluctuating']
+        variances = [v for v, _ in zip(iter(SymbolGen("variance")), lb_method.moments)] if params['fluctuating'] is True else params['fluctuating']
        correction = fluctuation_correction(lb_method, random_symbol(collision_rule.subexpressions, dim=lb_method.dim),
                                            variances)
+
        for i, corr in enumerate(correction):
            collision_rule.main_assignments[i] = Assignment(collision_rule.main_assignments[i].lhs,
                                                            collision_rule.main_assignments[i].rhs + corr)

+        temperature = sp.Symbol("fluctuating_temperature") if params["temperature"] is None else params["temperature"]
--- a/lbmpy/fluctuatinglb.py
+++ b/lbmpy/fluctuatinglb.py
 from lbmpy.moments import MOMENT_SYMBOLS
 from pystencils.simp.assignment_collection import SymbolGen
+from pystencils import Assignment
+
+
+def fluctuating_variance_equations(method, temperature=sp.Symbol("fluctuating_temperature"),
+                                   c_s_sq=sp.Symbol("c_s") ** 2,
+                                   variances=SymbolGen("variance")):
+    """Produces variance equations according to (3.54) in Schiller08"""
+    normalization_factors = abs(method.moment_matrix) * sp.Matrix(method.weights)
+    relaxation_rates_sqr = [r ** 2 for r in method.relaxation_rates]
+    density = method.zeroth_order_equilibrium_moment_symbol
+    mu = temperature * density / c_s_sq
+
+    return [Assignment(v, sp.sqrt(mu * norm * (1 - rr)))
+            for v, norm, rr in zip(iter(variances), normalization_factors, relaxation_rates_sqr)]