Opinion dynamics over complex networks: kinetic modeling and numerical methods

In this paper we consider the modeling of opinion dynamics over time dependent large scale networks. A kinetic description of the agents' distribution over the evolving network is considered which combines an opinion update based on binary interactions between agents with a dynamic creation and removal process of new connections. The number of connections of each agent influences the spreading of opinions in the network but also the way connections are created is influenced by the agents' opinion. The evolution of the network of connections is studied by showing that its asymptotic behavior is consistent both with Poisson distributions and truncated power-laws. In order to study the large time behavior of the opinion dynamics a mean field description is derived which allows to compute exact stationary solutions in some simplified situations. Numerical methods which are capable to describe correctly the large time behavior of the system are also introduced and discussed. Finally, several numerical examples showing the influence of the agents' number of connections in the opinion dynamics are reported.


Introduction
In recent years, the importance of large scale social networks has grown enormously and their study has raised lots of attentions, with the aim to understand how their structure and connections may influence the spread of opinions and ideas through human networks [1,6,15,16,26]. A major research topic is how to model the information exchange and, in particular, to understand and analyze the effects of interpersonal influence on processes such as opinion formation and creation and removal of new connections. The latter aspect is closely related to the construction of graph models for complex networks and has emerged as one of the most active research fields [2,7,8,21,27]. The empirical studies of technological and communication networks has been actively investigated thanks to a huge amount of data coming from the online platforms. From the theoretical point of view it is an unprecedented laboratory for testing the collective behavior of large populations of agents [9,30]. The need to handle with millions, and often billions, of vertices implied a considerable shift of interest to large-scale statistical properties of graphs.
In this context kinetic theory may play a major rule in designing effective models to characterize the statistical features of the opinion dynamics over such large collection of data. In particular, it can be used to analyze the so called stylized facts of the dynamics, like the asymptotic degree distribution of the connections in the network and the large time opinion behavior. To this aim, in this paper, we extend the kinetic model of opinion formation introduced in [29] to the case where each agent possesses a certain number of connections in the network. These connections evolve accordingly to a preferential attachment dynamics for the removal and creation of new connections. In this sense, the model here proposed fall in the general class of kinetic models for socio-economic problems where the dynamics of the model is influenced by additional characteristics of the agents, like personal conviction, leadership and knowledge [5,11,18,24,25].
In principle, the modeling proposed here is not limited to a particular kind of opinion dynamics and one can adapt other models developed in the literature [10,17,28] to evolve over the network by following the ideas presented in this paper. We mention here that recently opinion models have been considered in the context of optimal control in [4,5,6]. In a recent note [6] we faced the solution of an optimal control problem for a model of opinion dynamics described by a system of ordinary differential equations over an evolving network. More precisely we considered a network with a fixed number of vertices and edges which modifies its configuration of connections in time through a preferential attachment rewiring process.
A further contribution of the present manuscript is the development of numerical methods which are capable to describe correctly the large time behavior of the system. In particular we will focus on finite-difference schemes for the mean-field description of the opinion model over the network inspired by the well-known Chang-Cooper method [13,12,14,20]. We remark that, at variance with the standard Chang-Cooper method, the Fokker-Planck model considered here is nonlinear. Similar schemes for nonlinear Fokker-Planck equations have been previously introduced in [12,20].
The rest of the paper is organized as follows. In Section 2 we introduce the kinetic model and describe the evolution of the network of connections. The main properties of the network and the evolution of some macroscopic quantities, like the mean and the variance of the opinion over the network are discussed. Next in Section 3 we derive a Fokker-Planck model for the opinion dynamics under the classical quasi-invariant scaling. This permits to compute asymptotic stationary solutions of the opinion over the graph in some simplified situations. Section 4 is devoted to the construction of numerical methods for the above problems. Monte Carlo methods for the Boltzmann model and finite difference schemes for the Fokker-Planck model which are capable to describe correctly the steady states of the system are introduced. Finally in Section 5 several numerical examples illustrate our findings and show the behavior of the model. In separate Appendices we report proofs of related to the main properties of the network and to the positivity preservation property of the finite difference scheme.

The kinetic model
In this section we introduce a general mathematical model based on a kinetic description for the study of the opinion formation on a large evolving network.

Opinion dynamic
Let us consider a large system of agents interacting through a given network. We associate to each agent an opinion w, which varies continuously in a closed subset whose bounds denote two extreme and opposite opinions, and its number of connections c, as a discrete variable varying between 0 and the maximum number of connections allowed by the network. Note that this maximum number typically is a fixed value which is several orders of magnitude smaller then the size the network.
We are interested in the evolution of the density function where w ∈ I, I = [−1, 1] is the opinion variable, c ∈ C = {0, 1, 2, . . . , c max } is a discrete variable describing the number of connections and t ∈ R + denotes as usual the time variable. For each time t ≥ 0 we can compute the following marginal density which defines the evolution of the number of connections of the agents or equivalently the degree distribution of the network. In the sequel we assume that the total number of agents is conserved, namely cmax c=0 ρ(c, t) = 1.
The overall opinion distribution is defined likewise as the following marginal density function We express the evolution of the opinions by a binary interaction rule. From a microscopic point of view we suppose that the agents modify their opinion through binary interactions which depend on opinions and number of connections. If two agents with opinion and number of connections (w, c) and (w * , c * ) meet, their post-interaction opinion is given by where w, w * ∈ I = [−1, 1] denote the pre-interaction opinions and w , w * the opinion after the exchange of information between the two agents. Note that, in the present setting the compromise function P (·, ·; ·, ·) depends both on the opinions and on the number of connections of each agent. In (2.5) the nonnegative parameter η influences the compromise rate while ξ, ξ * are centered random variables with the same distribution Θ with finite variance ς 2 and taking values on a Borel set B ⊂ R. The function D(·, ·) ≥ 0 describes the local relevance of the diffusion for a given opinion and number of connections. We will consider by now a general interaction potential such that 0 ≤ P (w, w * , c, c * ) ≤ 1.
In absence of diffusion, ξ, ξ * ≡ 0, from (2.5) we have then the post-exchange distances between agents are still in the reference interval [−1, 1] if we consider η ∈ (0, 1/2) and 0 ≤ P (w, w * , c, c * ) ≤ 1. In agreement with [4,5,17,29] we can state the following result which derives the conditions on the noise term to ensure that the post-interaction opinions do not leave the reference interval.
Proposition 1. If we assume that 0 < P (w, w * ; c, c * ) ≤ 1 and then the binary interaction rule (2.5) preserves the bounds being the post interaction opinions w, w * contained in I = [−1, 1].
The evolution in time of the density function f (w, c, t) is described by the following integro-differential equation of Boltzmann-type where N [·] is an operator which is related to the evolution of the connections in the network and Q(·, ·) is the binary interaction operator defined as follows where ( w, w * ) are the pre-interaction opinions generated by the couple (w, w * ) after the interaction. The term J denotes the Jacobian of the transformation (w, w * ) → (w , w * ) and the kernels B, B define the binary interaction. Here and in the rest of the section, for notation simplicity, the explicit dependence from the time variable is omitted. We will consider interaction kernels of the following form where λ > 0 is a constant relaxation rate representing the interaction frequency.
In order to write the collision operator Q(·, ·) in weak form we consider a test function ψ(w) to get where the brackets < · > denotes the expectation with respect to the random variables ξ, ξ * . Equation (2.7) assumes the following weak form (2.11) An alternative form, obtained by symmetry is the following (2.12)

Evolution of the network
We introduced in the previous paragraph the operator N [·] characterizing the evolution of the agents in the discrete space of connections. This, of course, corresponds to the evolution of the underlaying network of connections between the agents. Here we will specify the details of the model considered in the present paper, inspired by [31]. The operator N [·] is defined through a combination of preferential attachment and uniform processes describing the evolution of the connections of the agents by removal and adding links in the network. These processes are strictly related to the generation of stationary scale-free distributions [8].
More precisely, for each c = 1, . . . , c max − 1 we define where γ = γ(t) is the mean density of connectivity defined as cρ(c, t), (2.14) α, β > 0 are attraction coefficients, and V r (f ; w) ≥ 0, V a (f ; w) ≥ 0 are characteristic rates of the removal and adding steps, respectively. The first term in (2.13) describes the net gain of f (w, c, t) due to the connection removal between agents whereas the second term represents the net gain due to the connection adding process. The factor 2 has been kept in evidence since connections are removed and created pairwise.
At the boundary we have the following equations which are derived from (2.13) taking into account the fact that, in the dynamics of the network, connections cannot be removed from agents with 0 connections and cannot be added to agents with c max connections.
Remark 1. If one defines the characteristic rates as and U a , U r constants, the dynamics in (2.13) correspond to a combination of a preferential attachment processes (α, β ≈ 0) and a uniform processes (α, β 1) for each agent with opinion w, with respect to the probability density of connections f (w, c, t)/g(w, t).
The evolution of the network of connections can be recovered taking ψ(w) = 1 in the master equation (2.11 Then for the collisional operator defined in (2.10) and the choice of N [·] in (2.13) the total number of agents is conserved. Let us take into account the evolution of the mean density of connectivity γ defined in (2.14). We can prove that, for each t ≥ 0 Therefore γ is not conserved in general. Asymptotically, conservation is recovered in the case β = 0, if the characteristic rates are given by (2.16) with U a = U r or are constants with V a = V r , and for a sufficiently fast decay of the density function f (w, c max , t).
In Appendix A.1-A.2 we report the explicit computations of the conservation of the total number of connections (2.19) and of the evolution of the mean density of connectivity (2.20).
In the particular case where V a and V r are constants independent of f and w then the operator N [·] is linear and will be denoted by L[·]. In this case, the evolution of the network of connections is independent from the opinion and we get the closed form and at the boundary (2.23) Note that the dynamics in (2.22) corresponds again to a combination of preferential attachment processes (α, β ≈ 0) and uniform processes (α, β 1) with respect to the probability density of connections ρ(c, t).
Concerning the large time behavior of the network of connections, in the linear case with V r = V a , β = 0 and now denoting by γ the asymptotic value of the density of connectivity, it is possible to prove that Proposition 2. For each c ∈ C the stationary solution to (2.21) or equivalently is given by (2.26) Detailed computations are given in Appendix A.3. Further approximations are possible in the cases α 1 and α ≈ 0. For big values of α the preferential attachment process described by the master equation (2.22) is destroyed and the network approaches to a random network, whose degree distribution is the Poisson distribution. In fact, in the limit α → +∞ we have (α + γ) c ≈ α(α + 1) · · · (α + c − 1) and In the second case, for γ ≥ 1 and small values of α, the distribution can be correctly approximated with a truncated power-law with unitary exponent

Evolution of the moments
In order to study the evolution of the mean opinion, i.e.  We obtain Of course, if the compromise function P (·, ·; ·, ·) is symmetric with respect to the pairs (w, w * ) and (c, c * ) we have conservation of the overall opinion on the network In addition, if the network operator is linear the evolution of the mean opinion obeys the same closed differential equation of the network of connections and therefore all the conclusions of the previous section hold true also for the mean opinion on the network. More generally we will consider compromise functions P (·, ·; ·, ·) with the following form where 0 ≤ H(·, ·) ≤ 1 represents the positive compromise propensity and 0 ≤ K(·, ·) ≤ 1 a function taking into account the influence of number connections in the opinion exchange process. Note that, in this case, even if we consider a symmetric compromise function H and a linear network operator we have 30) and the evolution of the mean opinion cannot be expressed in closed form due to the influence of the different connections that the agents possess. This is a fundamental difference compared to classical kinetic models of opinion [29].
In the case of the second moment of the opinion φ(w) = w 2 , if we assume a symmetric function P , by denoting which, in the case of a linear operator L[·] with P = 1 and in absence of noise D = 0, simplifies to (2.32) Equation (2.32) together with (2.21) and (2.28) form a closed system for the evaluation of the second order moment of the opinion.

Fokker-Planck modelling
In general it is difficult to obtain analytic results on the large time behavior of the opinion for the kinetic equation introduced in the previous section. A step toward the simplification of the analysis is the derivation of asymptotic states of the Boltzmann-type equation derived from a simplified Fokker-Planck-type models [24]. Here we recall briefly the approach usually referred to as the quasi-invariant opinion limit [5,11,29].

Derivation of the model
The idea is to rescale the interaction frequency λ, the interaction propensity η and the diffusion variance ς 2 at the same time, in order to to maintain asymptotically the memory of the microscopic interactions. Let su introduce the scaling parameter ε > 0 and consider the scaling The above scaling corresponds to the case where the interaction kernel concentrates on binary interactions producing very small changes in the agents' opinion but at the same time the number of interactions becomes very large. From a modelling point of view we require that the scaling (3.1) preserves the macroscopic properties of the kinetic system in the limit ε → 0, namely the evolution of the mean and the variance of the opinion derived in Section 2.3. The scaled equation (2.11) reads with scaled binary interactions given by where ξ is a centered random variable with variance εσ 2 . Since as ε → 0 we have w → w we can consider the Taylor expansion of ψ around w to get Inserting this expansion in the binary interaction term of (3.2) we obtain where R(ε) indicates the remainder, given by Therefore, the scaled binary interaction term reads By similar arguments of [29] it can be shown rigorously that R(ε) in (3.6) decays to zero in the limit ε → 0. Thus, as ε → 0 we recover Integrating backwards by parts equation (3.8) we obtain the following Fokker-Planck differential equation for the evolution of the opinions' distribution through the evolving network where

Stationary solutions
In this section we will show how in some cases it is possible to compute explicitly steady state solutions of the Fokker-Planck system (3.9). We restrict to linear operators L[·] and asymptotic solutions of the following form where ρ ∞ (c) is the steady state distribution of the connections (see Proposition 2) and From the definition of the linear operator L[·] we have L[ρ ∞ (c)] = 0, so stationary solutions of type (3.11) satisfy the following equation Under some simplifications we can analytically solve equation (3.13) as it has been shown in [5,29]. If we assume (2.29), i.e. P (w, can be written as follows and if we further assume that K(c, c * ) =K(c * ) is independent of c, and H(w, w * ) =H(w) independent of w * we have Therefore, on the support of ρ ∞ (c), stationary solutions can be derived from the following equation  which corresponds to the solution of the following ordinary differential equation where the constant C 0 is chosen such that the total mass of g ∞ is equal to one.
Some explicit examples are give below.
1. In the case H ≡ 1 and D(w) = 1 − w 2 the steady state solution g ∞ is given by 2. For H(w, w * ) = 1 − w 2 and D(w) = 1 − w 2 the steady state solution g ∞ is given by In Figure 2 as an example we report the stationary solution

Numerical methods
In this section we consider the development of numerical methods for the kinetic models studied in the previous sections. First we consider direct simulation Monte Carlo methods for the Boltzmann model (2.8) introduced in Section 2. Here the major difficulty is to consider a probabilistic interpretation of the dynamics induced by the network operator N [·], whereas the opinion interaction follows the standard binary sampling approach (see [24] for details). Next we consider the derivation of numerical schemes for the Fokker-Planck model (3.9) derived in Section 3. In particular we will focus on the construction of finite-difference methods which are capable to describe correctly the large time behavior of the model. To this aim we will consider a nonlinear version of the Chang-Cooper type discretization which has the nice feature of preserving the steady states and the non negativity of the numerical solution [14,20].

Direct simulation Monte Carlo
One of the most common approaches to solve Boltzmann-type equations is based Monte Carlo methods. Let us consider the initial value problem given by equation (2.7) with initial condition f (w, c, t = 0) = f 0 (w, c), the solution at time t n = n · ∆t n , n ≥ 1 is obtained as a composition of the solutions of the following problems: we first integrate the network term for all c ∈ C along the time interval [t n , t n+1 ] then we solve the interaction step f (w, c, 0) =f (w, c, t n ). (4. 2) The described process may be iterated in order to obtain the numerical solution of the initial equation at each time step. At variance with standard Monte Carlo methods for opinion dynamics, see for example [24], here we face the additional difficulty of the network evolution. In the sequel we describe the details of the Monte Carlo method for the network evolution in the simplified case N [·] = L[·].
Let f n = f (w, c, t n ) the empirical density function for the density of agents at time t n with opinion w ∈ [−1, 1] and connections c ∈ C. For a any given opinion w the solution of the transport step is given for each c > 0 and c < c max by with boundary conditions and temporal discretization such that .

end for
The collision step may be solved through binary interaction algorithm [3,24], where the basic idea is to solve the binary exchange of information described by (2.5), under the quasi-invariant opinion scaling (3.1). The time-discrete scheme reads where we have made explicit the dependence of Q(f, f ) on the frequency of interactions 1/ε and with Q + ε (f n , f n ) we denoted the gain part, namely it accounts the density of opinions gained at position w after the binary interaction (2.5). The collisional step (4.6) is a convex combinations of probability density under the time step constrain ∆t ≤ ε, which has to be coupled with (4.5). For further details on the algorithm we refer to [3,24].
In Figure 3 we show the two stationary states, already presented in Figure 2, computed through the Monte Carlo procedure just described, where we use N s = 2 × 10 4 samples to reconstruct the density and scaling parameter ε = 0.01 and ∆t = ε.

Chang-Cooper type numerical schemes
We consider the Fokker-Planck system (3.9) that we will rewrite in the form where and P[f ] is given by (3.10). The above equation is complemented with the initial data f (w, c, 0) = f 0 (w, c) and considered in the domain (w, c) ∈ I × C with zero flux boundary condition on w. Note that in the variable c the equation is in discrete form and therefore the discretization we will consider acts only on the opinion variable w.
Let us introduce a uniform grid w i = −1 + i∆w, i = 0, . . . , N with ∆w = 2/N , we denote by w i±1/2 = w i ± ∆w/2 and define Integrating equation (4.7) yields where F i [f ] is the flux function characterizing the numerical discretization.
We assume a flux function as a combination of upwind and centered discretization as in the classical Chang-Cooper flux where D i+1/2 = D(w i+1/2 , c) and D i+1/2 = D (w i+1/2 , c). The weights δ i+1/2 have to be chosen in such a way that a steady state solution is preserved. Moreover, as it is shown in Appendix B, this choice permits also to preserve nonnegativity of the numerical density.
Preservation of the steady states corresponds to assume that the numerical flux vanishes when f is at the steady state f ∞ . Imposing the numerical flux equal to zero from (4.10) we get Solving with respect to δ i+1/2 yields On the other hand the same computation directly on the flux (4.8) gives the differential equation which in general cannot be solved, except is some special cases as discussed in the previous section, due to the nonlinear term on the right hand side. A possible way to overcome this difficulty is to consider a quasi steady-state approximation as follows. We first integrate the previous equation in the cell [w i , w i+1 ] to get and then Next we can approximate the integral on the right hand side with a suitable quadrature formula. Because of singularities at the boundaries w = ±1 of the integrand function we can resort on open formula of Newton-Cotes type. For example, using the simple midpoint rule a second order approximation is obtained (4.14) Now by equating (4.14) and (4.11) we recover the following expression of the weight functions where (4.16) Note that here, at variance with the standard Chang-Cooper scheme [14], the weights depend on the solution itself as in [20]. Thus we have a nonlinear scheme which preserves the steady state with second order accuracy. In particular, by construction, the weight in (4.15) are nonnegative functions with values in [0, 1]. Higher order accuracy of the steady state can be recovered using a more general numerical flux given by and taking (4.18)

Numerical Results
In this section we perform several numerical test to validate our modeling and numerical setting. We focus on the case α < 1, since it represents the most relevant case in complex networks [2,31], for this range of the parameter we have emergence of power law distributions for the network's connectivity. Except for the first test case where we analyze the numerical convergence of the Boltzmann model in the quasi-invariant limit, in all the other test the opinion dynamics evolves according to (3.9). The compromise function P (c, c * ; w, w * ) and the local diffusion function D(w, c) will be specified in the various tests. The choice of parameters for the different tests is summarized in Table 1, where other parameters are introduced additional details will be reported.

Test #1
We first consider the simple test case where the opinion evolves independently of the network and we validate the Chang-Cooper type scheme comparing its convergence with respect to the Monte-Carlo methods. We simulate the dynamics with the linear compromise function P (w, w * ; c, c * ) = 1, and D(w, c) = 1 − w 2 , thus we can use the results (3.19), to compare the solutions obtained through the numerical scheme with the analytical one. The other parameters of the model are reported in Table 1 and we define the following initial data In Figure 4, on the left hand-side, we report the qualitative convergence of the Monte-Carlo methods, where we consider N s = 10 5 samples to reconstruct the opinion's density, g(w, t), on a grid of N = 80 points. The figure shows that for decreasing values of the scaling parameter ε = {0.5, 0.05, 0.005}, we have convergence to the reference solutions, (3.19) of the Fokker-Planck equation. On the right we report the convergence to the stationary solution of the connectivity distribution, (A.9), for α = 0.1 and V = 1 and with c max = 250. In this case we show two different qualitative behaviors for an increasing number of samples N s = {10 3 , 10 5 } and for sufficient large times.
In Figure 5, on the left-hand side we report the qualitative solution of the Chang-Cooper type scheme integrated with the explicit Euler method, on the right-hand side we depict the decay of the L 1 relative error to the reference solution, (3.19), i.e.
with g N representing the approximated solution of the numerical scheme. We test the scheme's convergence for different quadrature rules, and additionally we compared them with the error of the Monte-Carlo simulation. The right-hand plot shows how the Change Cooper scheme is able to reach high order of accuracy, for high order quadrature rules in (4.17). For the discretization we consider the following parameters: we account N = 80 and we define ∆w = 2/N and time step ∆t = (∆w) 2 /(4σ 2 ), for a final time of T = 10. The Binary Interaction algorithm is performed with ε = 0.0005 and with a number of sample N s = 10 5 .

Test #2
Next we consider a second validation test in the full case where the opinion and the network evolution are coupled, again with linear compromise function, P (w, w * ; c, c * ) = 1, and D(w, c) = 1 − w 2 . In this case we are able again to characterize the analytical solution of the model as the product of the two stationary solution for the opinion variable and the connectivity, i.e. f ∞ (w, c) = ρ ∞ (c)g ∞ (w). We define an initial data as follows, where and with coefficient c 0 = 20 and k 0 = 3/(20γ 3 0 ). Parameters are defined in Table 1, except for the characteristic rates, which will be defined in two different ways. We want to study the decay of the L 1 relative error with respect to the time, as depicted in Figure 5.
In the first case we consider constant characteristic rates, i.e. V = V a = V r , showing that for increasing values of V the convergence of the numerical scheme is faster. This is not surprising since for larger values of V the dynamics of the connectivity distribution relaxes faster towards the stationary state.
We report in Figure 6 the evolution of the density f (w, c, t), in the time frame [0, T ], with T = 20, where on the (z,c)-axis the distribution of the connections, ρ(c, t), is represented in order to better enlighten the convergence to the power-law like distribution.
In a second test we performed the same simulation, but with characteristic rates defined as in Remark 1, thus with U = U a = U r , β = 0 and γ f (t) = cmax c=0 cf (w, c, t). Simulations shows that in this case the same stationary solution are obtained.
We report in Figure 7, the decay of the errors for different values of the characteristic rates, in the two different cases, V = {10 3 , 10 4 , 10 5 } for the constant rate and  U = {10 3 , 10 4 , 10 5 } for the variable rates. In both cases we observe a faster convergence to the stationary solution for increasing values of the characteristic rates. Observe that the same order of accuracy of the mid-point rule in Figure 5 is recovered, on the other hand tiny differences in the decay are observed between the two cases. In Figure 8, we enlighten the different evolution of the transient solution at time t = 1, of the simulation in Figure  6. On the left we depict the solution with constant characteristic rates, on the right with variable characteristic rates, which shows that lower density in the opinion leads to faster spread on the connections.

Test #3
In this test we analyze the influence of the connections over the opinion dynamics, by considering a compromise function of the type where H(w, w * ) = 1 − w 2 and K(·, ·), which models the influence of the connectivity on the opinion evolution, is defined as follows for a, b > 0. This type of kernel assigns higher relevance into the opinion dynamics to higher connectivity, and low influence to low connectivity. The diffusivity is weighted by D(w, c) = 1 − w 2 . We perform a first computation with initial data where the parameters' choice is reported in the third line of Table 1, and a = b = 3 for the interaction function K(·, ·), (5.6). The evolution is performed through the Chang-Cooper type scheme with ∆w = 2/N , with N = 80. We study the evolution of the system in the time interval [0, T ], with T = 2.5.
In Figure 9 we report the result of the simulation. The initial configuration is is split in two parts, the majority concentrated around opinionw F = −1/2 and only a small On the left in the case of constant characteristic rate on the right variable characteristic rates defined as in (5.5). The right plot shows that for lower opinion's density the evolution along the connection is faster and slower where the opinions are more concentrated.
portion concentrated aroundw L = 3/4, from the upper-right and bottom-left figures we observe that, because of the anisotropy induced by K(c, c * ), the density with a low level of connectivity is immediately influenced by the small concentration of density around w L with a large level of connectivity; the bottom-right plot shows the final configuration. In Figure 10 we depict, on the left-hand side, the initial and final marginal density of the opinion, respectively g(w, 0) and g(w, T ), showing the change of the total opinion. On the right we enlighten the change of opinion plotting the evolution of the average opinion.

Test #4
In the last test case, we consider the Hegselmann-Krause model, [19], known also as bounded confidence model, where agents interact only with agents whose opinion lays within a certain range of confidence. Thus we define the following compromise function P (w, w * ; c, c * ) = χ {|w−w * |≤∆(c)} (w * ), (5.8) where ∆(c) is the confidence level and we assume that in general it depends on the number of connections. We define the initial data  Table 1 and D(w, c) = 1 − w 2 . The evolution is performed through the Chang-Cooper type scheme with ∆w = 2/N , with N = 80. We consider the evolution of the system in the time interval [0, T ], with T = 100. We study first a confidence level independent from the number of connections, therefore we set ∆(c) = ∆ = 0.25. In Figure 11 the evolution of the initial data (5.9) shows the classical behavior of Hegselmann-Krause model, where opinions' clusters emerge.
Next, we perform a second computation where the confidence bound depends on the number of connections as follows This choice reflects a behavior where agents with higher number of connections are prone to larger level of confidence. We report in Figure 12 the evolution of (5.9), where ∆(c) creates an heterogeneous emergence of clusters with respect to the connectivity level: for higher level of connectivity consensus is reached, since the bounded confidence level is larger, instead for lower levels of connectivity multiple clusters appears, up to the limiting

Conclusions
The construction of kinetic models and numerical methods for the spreading of opinions over time dependent large scale networks has been considered. First we have introduced a Boltzmann model for the opinion interactions based on a preferential attachment process for the creation of new connections between agents. If the preferential attachment is independent from the agents' opinion the large time behavior of the network can be described analytically and originates both Poisson type distributions as well as truncated power laws. Next we derived the corresponding mean-field approximation which permits to have a deeper understanding of the asymptotic behavior of the opinion dynamics and to compute analytically stationary states in simplified situations. Robust numerical methods, based on stochastic as well as deterministic techniques have been introduced and their property discussed. The results, for various test cases, show the validity of the present approach. Several extensions of the present approach are possible. First one may consider the case where the number of agents in the network is not conserved, as it happens in a real social network. Moreover, control problems with the aim to force consensus over the network may be introduced. From (2.13) we have Using the boundary conditions (2.15) we have the desired property. As a consequence we obtain the conservation of the total number of agents A.2 Mean density of connectivity Next we consider the evolution of the mean density of connectivity γ(t). We prove that In fact, thanks to (2.13), in the internal points we have We observe that the first sum in (A.4) is equal to Similarly the second sum in (A.4) is equal to (A.6) Using the boundary condition for c = c max (since the one at c = 0 does not play any role here) we have which together with the above computations yields (A.3).
As a consequence we have

A.3 Asymptotic behavior
In the following we compute the explicit stationary solution ρ ∞ (c) for the evolution of ρ(c, t) in the linear case with V a = V r , β = 0 and assuming Note that in the sequel, for notation simplicity, we denote by γ = γ ∞ the asymptotic stationary value reached by the mean density of connectivity.
Proposition 3. For each c ∈ C the stationary solution to (2.21) or equivalently is given by (A.11) Proof. Let us show (A.10) by induction. First, from the boundary condition (2.23) at c = 0 we immediately have Now let us assume that (A.10) holds true for c, we want to prove that By direct inspection one verifies that also the boundary condition (2.23) at c = c max is verified.

B Properties of the implicit-explicit scheme
Let us consider the following implicit-explicit discretization of (3.9) where f n i = f n i (c), endowed with a positive initial condition f 0 i (c) = f i (c, 0). The main motivation for the time discretization above is related to the severe stability constraints of an explicit scheme applied to the network operator which would require the time step to be O(1/c max ) where c max 1.

B.1 Positivity
In order to study the nonnegativity property of scheme (B.1) it is convenient to rewrite it as a sequence of two steps The first step involves the Chang-Cooper type scheme and reads where B n i+1/2 , C i+1/2 are given by In fact, setting x = B n i+1/2 ∆w/C i+1/2 , y = B n i−1/2 ∆w/C i−1/2 the two inequalities are equivalent to show that ∀ x, y ∈ R Typically when σ 2 is large this will originate a parabolic stability condition that requires ∆t = O(∆w 2 ). This can be avoided taking the diffusive part implicitly, however, since we were mostly interested in the case of small values of σ 2 we will not pursue this direction here.
Next, we consider the second step

B.2 Conservations and stability
Let us consider the conservation properties of the scheme with respect to the variable w. Let us observe that from scheme (B.1) we get