SCALING AND NUMERICAL METHOD

Theorem 4.6.1. In R², letp > 1 and assume that n₀ ∈ L¹₊(R², (1 + |x|²)dx). There exists a constantα such that when the initial data satisfies m₀(0) < _{p A S κ}^4α

c, there exists a weak solution to (4.43) inL^p(R², dx) at all times.

Proof. The estimates are derived based on the Sobolev inequalities following the argu-ment by Jager and Luckhaus [123].

Multiplying (4.43) by n^p−1and integrating, one can get 1

To estimate the integral with power p + 1, the following Gagliardo-Nirenberg-Sobolev inequality is employed

R²n^pdx decays in time. From this a priori estimate, one may conclude the existence as done in [127].

4.7 Scaling and numerical method

The present section and the following section deal with the numerical study of the model (4.8). This section presents the dimensionless form of the model and the finite ele-ment method employed to discretize the differential equations in (4.8). All computations presented are realized by means of the FreeFem++ code presented in chapter 2.

4.7.1 Scaling

The dimensionless variables n_g⁰, n_d⁰, c⁰, s_n⁰, s_e⁰, x⁰, and t⁰are defined by n_g⁰ = n^refn_g, n_d⁰ = n^refn_d, c⁰ = c^refc, s_n⁰ = s^ref_n s_n,

s_e⁰ = s^ref_e s_e, x⁰ = ¹_{`</sub>x, t<sup>0</sup> = <sup>D</sup><sub>`}^m2 t. (4.53)

4.7. SCALING AND NUMERICAL METHOD

Thus, introducing (4.53) in (4.8) and removing the prime of the dimensionless variables, the dimensionless system given below can be derived

∂n_g

∂t − ∇²n_g+ Π_S∇. (n_g∇c) = −Π_AΦn_g+ Π_krn_d, t > 0, x ∈ Ω

∂c

∂t − Π_Dc∇²c = Π_κcΦn_g− Π_δcc,

∂n_d

∂t − ∇²n_d= Π_AΦn_g− Π_krn_d,

∂s_n

∂t − Π_Dsn∇²s_n= Π_κnΦn_g− Π_δsns_n,

∂s_e

∂t − Π_Dse∇²s_e = Π_κeΦn_g− Π_δses_e.

(4.54)

where

ΠS= S c^ref

D_m , ΠA = ^`_D^2A

m, Π_kr = `²k_r D_m , Π_Dc = D_c

D_m, Π_κc = ^κc_D^A`2nref

mc^ref , Π_δc = `²δ_c D_m, Π_Dsn = D_sn

D_m, Π_κn = ^κnA`2nref

Dms^refn

, Π_δsn = `²δ_sn D_m , Π_Dse = D_se

D_m, Π_κe = ^κeA`2nref

Dms^refe , Π_δse = `²δ_se D_m . The reference quantities are chosen as

n^ref = initial concentration of granulated mastocytes,

c^ref = concentration of the stored chemoattractant before degranulation, s^ref_n = concentration of the stored nerve stimulant before degranulation, s^ref_e = concentration of the stored endocrine stimulant before degranulation.

4.7.2 Finite element method

For the sake of clarity, let us simplify the notation by changing the previous dimen-sionless parameters Π_iinto i. Then, the chemotaxis system is solved in the domain Ω in a

4.7. SCALING AND NUMERICAL METHOD

The above system is endowed with the prescribed initial conditions in (4.9) and the homogeneous Neumann boundary condition (4.10) or the zero total flux boundary condi-tions (4.11).

4.7.2.1 Weak formulation

Let the solution (ng, n_d, c, s_n, s_e) of (4.55) be regular enough, for example in (H²(Ω))⁵. The equations in (4.55) are multiplied by a test function w_iin (H¹(Ω))⁵(i = g, d, c, n, e) and are integrated over Ω. Applying the Green formula and taking into account the homo-geneous Neumann boundary condition (4.10) yield the weak formulation of the problem in the subspace H¹(Ω) that reads as follows. Find (n_g, n_d, c, s_n, s_e) ∈ (H¹(Ω))⁵ such as

The term n_g∇c is not easy to treat because of its nonlinear nature when considering a monolithic method. It is convenient to decouple the system using the previously com-puted solutions. For the temporal discretization scheme, the θ - scheme given in (??) is employed.

4.7. SCALING AND NUMERICAL METHOD

A semi-discretization in time gives the following scheme.

1. Solve the semi-linear scheme in c^m+1using n^m_g c^m+1− c^m

∆t − D_c∇²c^m+1 = κ_cΦn^m_g − δ_cc^m+1. (4.57) The solution c^m+1 serves to compute the mastocyte density.

2. Solve the semi-linear scheme in n^m+1_d using n^m_g n^m+1_d − nⁿ_d

∆t − ∇²n^m+1_d = AΦn^m_g − k_rn^m+1_d . (4.58) The solution n^m+1_d serves to compute the mastocyte density.

3. Solve the semi-linear scheme in n^m+1_g using c^m+1, n^m+1_d n^m+1_g − n^m_g

∆t − ∇²n^m+1_g + S∇ · n^m+1_g ∇c^m+1
= −AΦn^m+1_g + k_rn^m+1_d . (4.59) 4. Solve the semi-linear scheme in s^m+1_n using n^m+1_g

s^m+1_n − s^m_n

∆t − D_sn∇²s^m+1_n = κ_nΦn^m+1_g − δ_sns^m+1_n . (4.60) 5. Solve the semi-linear scheme in s^m+1_e using n^m+1_g

s^m+1_e − s^m_e

∆t − D_se∇²s^m+1_e = κ_eΦn^m+1_g − δ_ses^m+1_e . (4.61) 4.7.2.3 Full discretization

The full discretization of the problem (4.56) corresponds to the finite element space discretization of the semi-discretized scheme represented by (4.57)-(4.61). The scheme for the model (4.8)-(4.9)-(4.10) is summarized as follows

1. Find c^m+1 ∈ V_c^h such that

4.7. SCALING AND NUMERICAL METHOD

Remark 3. To complete the numerical method, one needs to choose the finite element spaces. The choice of the classicalP¹ finite element space forn^m+1_g , n^m+1_d , s^m+1_n , and s^m+1_e is reasonable. However, forcⁿ⁺¹ the classicalP² finite element space is required.

Indeed, ifc^m+1is piecewise P2, then the chemotactic vector fieldS∇c^m+1is piecewise P1. The scheme can be applied to the models (4.36)-(4.9)-(4.10) and (4.43)-(4.9)-(4.10) by skipping the steps corresponding to the omitted equations.

One of the advantages of the artificial decoupling is the quasi-linearity of the equa-tions. This reduces the computational cost and thus the efficiency of the method is im-proved. Nevertheless, the quality of the solutions depends on the chemotactic sensitiv-ity S, the activation rate A, and the release coefficient κ_c. To enhance the quality of the computation (and to reduce the computation cost), mesh adaptation (a subroutine of FreeFem++) fits the initial condition, i.e., a given cell distribution within the domain of interest, as the solution evolves locally.

4.7.2.4 Code validation

Verification of the traditional Keller-Segel equation To quantify the convergence be-havior with respect to the mesh size, several meshes are considered on the domain Ω and the errors are evaluated on the density n_g and the concentration c with respect to the L²(Ω) and H¹(Ω) norms. The related convergence are plotted in figures 4.2, 4.3, and 4.4 in logarithmic scales. For the P1-P2 finite elements, the slope for n_g in the L²(Ω) norm is evaluated to be 2.50, in the L^∞(Ω) norm evaluated to be 2.11 and in the H¹(Ω) norm

4.7. SCALING AND NUMERICAL METHOD

evaluated to be 1.01. The slope for c in the L²(Ω) norm is evaluated to be 2.11, in the L^∞(Ω) norm evaluated to be 2.10 and in the H¹(Ω) norm evaluated to be 2.08. For the P2-P2 finite elements, the slope for n_g in the L²(Ω) norm is evaluated to be 2.57, in the L^∞(Ω) norm evaluated to be 2.70 and in the H¹(Ω) norm evaluated to be 2.37. The slope in c in the L²(Ω) norm is evaluated to be 2.06, in the L^∞(Ω) norm evaluated to be 2.07 and in the H¹(Ω) norm evaluated to be 2.07. The rate of convergence for n_gin the H¹(Ω) norm is considerably improved. For the P2-P3 finite elements, the slope for n_g and c are close to the previous one. Thus for the P2-P3 finite elements, the rates of convergence for n_gand c are not improved. From the computed convergence rates the choice of P2-P2

finite elements drastically improves the quality of the solution. Therefore, from now on, numerical simulations will be carried out using P2-P2 finite elements

1e-06 1e-05 0.0001 0.001 0.01 0.1 1 10

0.001 0.01 0.1 1

∆x

L^∞norm (ng) L²norm (ng) H¹norm (ng) L^∞norm (c) L²norm (c) H¹norm (c)

Figure 4.2: Convergence rates for the P1-P2 finite elements

Sensitivity to the chemotactic sensitivity The next issue is to check the sensitivity of the computation to the chemotactic sensitivity parameter S. Numerical simulations are carried out with a fixed mesh where the chemotactic sensitivity parameter S varies in the interval [1, 100]. The error curves are then plotted against S for different meshes. Fig-ures 4.5, 4.6, and 4.7 show the curves for mesh sizes ∆x = 1/16, ∆x = 1/64, and

∆x = 1/256, respectively. The important observation is that the chemotactic sensitivity parameter affects strongly the results of n_g. However, the chemotactic sensitivity param-eter seems to affect very slightly the results of c. This shows the necessity to refine the mesh to improve the quality of the results. Mesh adaptation can be used to refine the mesh locally according to the weight of the chemotaxis term S∇ · (n_g∇c).

4.7. SCALING AND NUMERICAL METHOD

1e-06 1e-05 0.0001 0.001 0.01 0.1 1

0.001 0.01 0.1 1

∆x

L^∞norm (ng) L²norm (ng) H¹norm (ng) L^∞norm (c) L²norm (c) H¹norm (c)

Figure 4.3: Convergence rates for the P2-P2 finite elements

1e-06 1e-05 0.0001 0.001 0.01 0.1 1

0.001 0.01 0.1 1

∆x

L^∞norm (ng) L²norm (ng) H¹norm (ng) L^∞norm (c) L²norm (c) H¹norm (c)

Figure 4.4: Convergence rates for the P2-P3 finite elements

4.7. SCALING AND NUMERICAL METHOD

1e-06 1e-05 0.0001 0.001 0.01 0.1 1 10

0.1 1 10 100 1000 10000

S L^∞norm (ng)

L²norm (ng) H¹norm (ng) L^∞norm (c) L²norm (c) H¹norm (c)

Figure 4.5: Influence of the chemotactic sensitivity parameter S with a mesh size ∆x = 1/16.

1e-06 1e-05 0.0001 0.001 0.01 0.1 1 10

0.1 1 10 100 1000 10000

S L^∞norm (ng)

L²norm (ng) H¹norm (ng) L^∞norm (c) L²norm (c) H¹norm (c)

Figure 4.6: Influence of the chemotactic sensitivity parameter S with a mesh size ∆x = 1/64.

4.8. COMPUTATIONAL MODEL

1e-06 1e-05 0.0001 0.001 0.01 0.1 1 10

0.1 1 10 100 1000 10000

S L^∞norm (ng)

L²norm (ng) H¹norm (ng) L^∞norm (c) L²norm (c) H¹norm (c)

Figure 4.7: Influence of the chemotactic sensitivity parameter S with a mesh size ∆x = 1/256.

4.8 Computational model

4.8.1 Acupoints

Acupoints feature a higher mastocyte density than its neighboring regions. The initial mastocyte distribution can be modeled by one Gaussian distribution given below

n(x, y, t = 0) = a_nexp − (x − x_n0)²

2σ²_nx + (y − y_n0)² 2σ_n²_y

!!

+ b_n, (x, y) ∈ R². (4.67)

In the above, a_n is the amplitude of the Gaussian distribution, b_nthe minimal mastocyte density in the tissue, σ_nxthe dispersion in the x-direction, σ_ny the dispersion in the y-direction, and (x_n0, y_n0) the coordinates of maximum mastocyte density.

In the full space R², if b_n = 0, σ_nx = σ_ny = σ, and x_n0 = y_n0 = 0, the number of mastocyte is given by

m_|

R2 = Z _−∞

−∞

Z _−∞

−∞

n(x, y, t = 0)dxdy = 2πa_nσ². (4.68)

4.8. COMPUTATIONAL MODEL

On the disc centered at (0, 0) with radius R, the number of mastocyte is given by m_|D(0,R) =

By choosing accordingly the parameter a_n, σ_nx and σ_ny, the quantity of mastocytes and their dispersion can be controlled (see figure 4.8.1 and table 4.2) in the initial distri-bution. Acupointa are represented by a large quantity and the concentrated distribution of mastocytes, whereas non-acupoints are characterized by a dispersed distribution of mas-tocytes (see figure 4.10).

Table 4.2: Repartition of the mass in a Gaussian distribution in the full domain R and R²

1D 2D

The stress function Φ is smooth and compactly supported as defined in (4.12). In the numerical simulations, the so-called bump function defined below is employed

Φ_B(x, y) =

In the above, ` is the distance range of the applied stress. α and β are the two positive parameters that control the amplitude and the shape of the function Φ (see figure 4.9)

4.8. COMPUTATIONAL MODEL

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 4.8: Gaussian spatial distribution of mastocytes in a bounded domain Ω. If b_n= 0 and σ_nx = σ_ny = σ = ^R₂, then approximately 86.5% of the total initial mastocytes is in the disc D((x_n0, y_n0), 2σ).

-1 -0.5

0 0.5

1 -1 -0.5

0 0.5

1 0.10

0.20.3 0.40.5 0.60.7 0.80.91

α = 2, β = 1 α = 1, β = 0.1 α = 1e5, β = 10

Figure 4.9: Stress function Φ defined as the bump function (4.70)

在文檔中 Modeling and simulation of transport during acupuncture (頁 114-125)