6. Results and Discussion
6.7. PDE Solution
7 // d i f f u s i o n
8 a :A <−E−> b : A, ( a < b ) : :> {
9 double r x = D X ∗ ( a [ x ] − b [ x ] ) ;
10 a [ x ] : ’= −r x ;
11 b [ x ] : ’= +r x ;
12 double r y = D Y ∗ ( a [ y ] − b [ y ] ) ;
13 a [ y ] : ’= −r y ;
14 b [ y ] : ’= +r y ;
15 }
Instead of a regular grid, one could also investigate how Turing patterns emerge on an arbitrary network. This idea has been pursued in [NM10].
6.7. PDE Solution
Let us recall the example of heat conduction in a metal stick from section 2.12.5. The PDEut = uxx models heat transport with an initial condition u(x,0) = sin(x·π/2) and boundary conditionsu(0, t) = 0 and ux(L, t) = 0.
The method of lines is applied to the stick with a spatial discretization into 21 points.
Each point is represented in XL by a moduleAdefined as
module A(double u ) extends Box ==> s e t C o l o r ( u , 0 , 0 ) ;
The attributeustores the current value of the temperature for this segment, which is visualized by a shade of red set using an instantiation rule. For neighbouring segments the corresponding nodes are connected by successor edges.
Using finite differences and taking into account the boundary conditions, the following set of ODEs must be integrated:
du0
dt = 0, dui
dt = ui−1−2ui+ui+1
∆x2 , 0< i < N duN
dt = 2·(uN−1−uN)
∆x2 ,
where N is the number of points, ∆x = L/N the spacing between the points, and ui
the temperature at point numberi. A central finite difference approximation of second order has been used for all points, except those at the borders. Foru0the first boundary
6.7. PDE Solution
Figure 6.12.: Reaction-diffusion system in two dimensions. Amount of morphogen in each cell is indicated by RGB colors (red component for morphogen X, green component for morphogenY).
Listing 6.4: XL rules used to apply method of lines to heat conduction.
1 a :A : :> {
2 A l = f i r s t ( (∗ a < A ∗) ) ;
3 A r = f i r s t ( (∗ a > A ∗) ) ;
4 i f ( l == n u l l) {
5 // no change , k e e p a t z e r o
6 } e l s e i f ( r == n u l l) {
7 a [ u ] : ’= 2 ∗ ( l [ u ] − a [ u ] ) / ∆x2;
8 } e l s e {
9 a [ u ] : ’= ( l [ u ] − 2∗a [ u ] + r [ u ] ) / ∆x2;
10 }
11 }
condition requires the value to remain zero, which is also true for its time-derivative.
For the other enduN of the stick we need the non-existent value uN+1, which can be derived from the finite difference approximation of the second boundary condition
∂uN
∂x ≈ uN+1−uN−1
2∆x = 0
asuN+1 =uN−1.
Listing 6.4 shows the equivalent formulation of the ODEs in XL. For each nodeathe left and right neighbour is searched in the graph. If the node was on one of the ends of the metal stick, such a neighbour is missing and this is handled accordingly.
Numerical integration of the model was done with an implicit Adams-Moulton inte-grator. The resulting time series is shown in figure 6.13. Note that explicit methods might have difficulties according to the CFL condition, especially if a finer discretization is used.
t = 0.0 : t = 0.1 : t = 0.2 : t = 0.3 : t = 0.4 : t = 0.5 : t = 0.6 : t = 0.7 : t = 0.8 : t = 0.9 : t = 1.0 :
Figure 6.13.: Method of lines solution for a metal stick.
In the modelling of trees, which currently is one of the main application fields of PDEs on Tree
Structures
6.7. PDE Solution
GroIMP/XL, segments of the trunk or the branches are often represented by cylinders as nodes in the graph. When accurate information about the spatial distribution of chemical substrates within the segments is needed, these segments must be discretized.
This raises the question whether there is an automatic way to achieve this.
If we take a look at figure 6.14a we can see that the structure, as it is usually mod-elled, causes segments to overlap spatially. This causes problems if the discretization is one-dimensional along the main axis of each segment (see figure 6.14b), since in the overlapping region it must be decided for each point in space to which segment it should be assigned. A three-dimensional discretization would perhaps solve this issue, but now we have way more ODEs to integrate. Also it is not clear how the inner workings of the transport mechanism should be accounted for (i.e. transport pipes within branches).
(a) overlapping segments (b) discretization
Figure 6.14.: Original tree structure (a) and a one-dimensional discretization (b) of it.
For these and other reasons, an automatic discretization of 3D structures is not feasi-ble, thus the user should be responsible to provide a proper one. However, there should be some way to support the user in doing this. Since the topology is already available, the user can attach to it the discretized structure that is created by ordinary XL rules.
This also gives freedom in interpretation of the topology, since this is fully under the user’s control.
Each segment of the tree structure (cf. figure 6.15a) can be decomposed like it was done for the metal stick before, by creating a chain of nodes connected by unidirectional edges.
These chains can then be linked to the original segment using some user-defined edge type, thus allowing to handle each segment separately in a first pass (see figure 6.15b.
Then, in a second pass, these chains can be linked together according to the topology
F
Figure 6.15.: Tree structure with attached spatial discretization.
(see figure 6.15c).
The program shown in listing 6.5 demonstrates how this idea can be implemented in XL. In init the structure from figure 6.15a is created. Then, in discretize, the first rule attaches a chain of user-defined nodes of typeAto each segmentF. The number of A nodes created could optionally depend on the length of theFsegment to obtain a more regular discretization, but for sake of simplicity a fixed number of three As was used here. Calling derive makes the modifications to the graph visible, so that the second rule can finally connect the chains according to the topology.
Usually, the structure created to represent the tree also contains other commands like rotations, translations, etc. and branches are connected by branch edges and not successor ones. So if the structure that ought to be discretized is created by this rule
Axiom ==> F F [RU( 3 5 ) F ( 1 . 5 ) ] RU(−30) F ;
we have to slightly modify our discretization function to take care of this. The second rule would then be replaced by
(∗ a1 :A −E−> f 1 : F (−( s u c c e s s o r|branch)−>)+ : ( f 2 : F) −E−> a2 :A ∗)
==>> a1 a2 ;
to ensure thatsuccessor and branchedges are properly skipped.