Carry out a few single indivual steps of diffusion,¶

and directly verify that the values satisfy the diffusion equation¶

In this "PART 1", we'll be looking at a tiny system with just 10 bins, to easily inspect the numbers directly. We'll use concentrations initialized to an upward-shifted sine wave with 1 cycle over the length system

TAGS : "diffusion 1D", "under-the-hood"¶

In [1]:

LAST_REVISED = "May 3, 2025"
LIFE123_VERSION = "1.0.0rc3"       # Library version this experiment is based on

In [2]:

#import set_path              # Using MyBinder?  Uncomment this before running the next cell!

In [3]:

#import sys
#sys.path.append("C:/some_path/my_env_or_install")   # CHANGE to the folder containing your venv or libraries installation!
# NOTE: If any of the imports below can't find a module, uncomment the lines above, or try:  import set_path   

from life123 import BioSim1D, ChemData, CollectionArray, Numerical, check_version

import numpy as np

In [4]:

check_version(LIFE123_VERSION)

OK

In [ ]:

In [5]:

# We'll be considering just 1 chemical species, "A"
diffusion_rate = 10.

chem_data = ChemData(diffusion_rates=[diffusion_rate], names=["A"])

In [6]:

# Initialize the system with just a few bins
bio = BioSim1D(n_bins=10, chem_data=chem_data)

In [7]:

# Initialize the concentrations to a sine wave with 1 cycle over the system
# (with a bias to always keep it > 0)
bio.inject_sine_conc(chem_label="A", number_cycles=1, amplitude=10, bias=50)

In [8]:

bio.show_system_snapshot()

SYSTEM SNAPSHOT at time 0:

Out[8]:

	A
0	50.000000
1	55.877853
2	59.510565
3	59.510565
4	55.877853
5	50.000000
6	44.122147
7	40.489435
8	40.489435
9	44.122147

In [9]:

# Visualize the initial system state
bio.visualize_system(title_prefix="Initial System State (for the tiny system)")

In [10]:

# Show as heatmap
bio.system_heatmaps(title_prefix="Initial System State (for the tiny system)")

In [ ]:

Now do 4 rounds of single-step diffusion, to collect the system state at a total of 5 time points:¶

t0 (the initial state), plus t1, t2, t3 and t4¶

In [11]:

# All the system states will get collected in this object
history = CollectionArray()

In [12]:

# Store the initial state
arr = bio.lookup_species(chem_index=0, copy=True)
history.store(par=bio.system_time, data_snapshot=arr, caption=f"State at time {bio.system_time}")

In [13]:

# Take a look at what got stored so far (a matrix whose only row is the initial state)
history.get_array()

Out[13]:

array([[50.        , 55.87785252, 59.51056516, 59.51056516, 55.87785252,
        50.        , 44.12214748, 40.48943484, 40.48943484, 44.12214748]])

In [14]:

# Additional parameters of the simulation run (the diffusion rate got set earlier)
delta_t = 0.01
delta_x = 2       # Note that the number of bins also define the fraction of the sine wave cycle in each bin
algorithm = None  # "Explicit, with 3+1 stencil"

In [15]:

# Do the 4 rounds of single-step diffusion; show the system state after each step, and accumulate all data
# in the history object
for _ in range(4):
    bio.diffuse(time_step=delta_t, n_steps=1, delta_x=delta_x , algorithm=algorithm)
    bio.describe_state(concise=True)

    arr = bio.lookup_species(chem_index=0, copy=True)
    history.store(par=bio.system_time, data_snapshot=arr, caption=f"State at time {bio.system_time}")

SYSTEM STATE at Time t = 0.01:
[[50.14694631 55.82172403 59.41974735 59.41974735 55.82172403 50.
  44.17827597 40.58025265 40.58025265 44.03132966]]
SYSTEM STATE at Time t = 0.02:
[[50.28881576 55.76980517 59.32979676 59.32979676 55.76613151 50.
  44.23386849 40.67020324 40.66652958 43.94505274]]
SYSTEM STATE at Time t = 0.03:
[[50.42584049 55.72178022 59.24079697 59.24070513 55.71106985 50.
  44.28893015 40.75920303 40.7485845  43.86308966]]
SYSTEM STATE at Time t = 0.04:
[[50.55823898 55.67735715 59.15281926 59.15246655 55.65653399 50.
  44.34346372 40.84718074 40.82671259 43.78522703]]

In [16]:

# Now, let's examine the data collected at the 5 time points
all_history = history.get_array()
all_history

Out[16]:

array([[50.        , 55.87785252, 59.51056516, 59.51056516, 55.87785252,
        50.        , 44.12214748, 40.48943484, 40.48943484, 44.12214748],
       [50.14694631, 55.82172403, 59.41974735, 59.41974735, 55.82172403,
        50.        , 44.17827597, 40.58025265, 40.58025265, 44.03132966],
       [50.28881576, 55.76980517, 59.32979676, 59.32979676, 55.76613151,
        50.        , 44.23386849, 40.67020324, 40.66652958, 43.94505274],
       [50.42584049, 55.72178022, 59.24079697, 59.24070513, 55.71106985,
        50.        , 44.28893015, 40.75920303, 40.7485845 , 43.86308966],
       [50.55823898, 55.67735715, 59.15281926, 59.15246655, 55.65653399,
        50.        , 44.34346372, 40.84718074, 40.82671259, 43.78522703]])

Each row in the above matrix is a state snapshot at a different time:¶

first row is the initial state at time 0; the successive rows are at t1, t2, t3, t4

In [17]:

# Let's consider the state, i.e. the concentration across the bins, at the midpoint in time (t2)
f_at_t2 = all_history[2]
f_at_t2

Out[17]:

array([50.28881576, 55.76980517, 59.32979676, 59.32979676, 55.76613151,
       50.        , 44.23386849, 40.67020324, 40.66652958, 43.94505274])

If one compares the above state (at t2) with the initial state (the top row in the matrix), one can see, for example, that the leftmost bin's concentration is increasing ("pulled up" by its neighbor to the right): 50. has now become 50.28881576

The rightmost bin's concentration is decreasing ("pulled down" by its neighbor to the left): 44.12214748 has now become 43.94505274

The diffusion equation states that the partial derivative of the concentration values with respect to time must equal the (diffusion rate) times (the 2nd partial derivative with respect to space).¶

Let's see if that is the case for the values at time t2! We are picking t2 because we need at least a value at the earlier time, and a value at the later time

A. The 2nd partial derivative with respect to space¶

In [18]:

# A simple-minded way of computing the 2nd spacial derivative is to use the Numpy gradient function TWICE across the x value
# (that function approximates the derivative using differences)
gradient_x_at_t2 = np.gradient(f_at_t2, delta_x)   # Start with taking the first derivative
gradient_x_at_t2    # This will be the partial derivative of the concentration with respect to x, at time t2

Out[18]:

array([ 2.74049471,  2.26024525,  0.8899979 , -0.89091631, -2.33244919,
       -2.88306575, -2.33244919, -0.89183473,  0.81871237,  1.63926158])

For example, the 2nd entry in the above array of estimated derivatives, can be manually checked from the 1st and 3rd value in the function f_at_t2(x), as follows:

In [19]:

(59.32979676 - 50.28881576) / (2*delta_x)    # This way of numerically estimating derivatives is called "Central Differences"

Out[19]:

2.2602452500000005

Now take the derivative again, with the Numpy gradient function, to arrive at a coarse estimate of the 2nd derivative with respect to x:¶

In [20]:

second_gradient_x_at_t2 = np.gradient(gradient_x_at_t2, delta_x)
second_gradient_x_at_t2

Out[20]:

array([-0.24012473, -0.4626242 , -0.78779039, -0.80561177, -0.49803736,
        0.        ,  0.49780776,  0.78779039,  0.63277408,  0.4102746 ])

Note how the 2nd derivative is 0 at bin 5 (bins are numbered 0 thru 9): if you look at the earlier sine plot, x5 is the inflection point.

B. The partial derivative with respect to time¶

Now, let's look at how concentrations change with time. Let's first revisit the full history:

In [21]:

all_history

Out[21]:

array([[50.        , 55.87785252, 59.51056516, 59.51056516, 55.87785252,
        50.        , 44.12214748, 40.48943484, 40.48943484, 44.12214748],
       [50.14694631, 55.82172403, 59.41974735, 59.41974735, 55.82172403,
        50.        , 44.17827597, 40.58025265, 40.58025265, 44.03132966],
       [50.28881576, 55.76980517, 59.32979676, 59.32979676, 55.76613151,
        50.        , 44.23386849, 40.67020324, 40.66652958, 43.94505274],
       [50.42584049, 55.72178022, 59.24079697, 59.24070513, 55.71106985,
        50.        , 44.28893015, 40.75920303, 40.7485845 , 43.86308966],
       [50.55823898, 55.67735715, 59.15281926, 59.15246655, 55.65653399,
        50.        , 44.34346372, 40.84718074, 40.82671259, 43.78522703]])

For simplicity, let's start by just inspecting how the values change over time at the 3rd bin from the left,¶

i.e. the 3rd column (index 2 because counting starts at 0) of the above matrix:

In [22]:

f_of_t_at_x2 = all_history[ : , 2]  # The column with index 2
f_of_t_at_x2

Out[22]:

array([59.51056516, 59.41974735, 59.32979676, 59.24079697, 59.15281926])

The above a function of time¶

(we took an entry from each row - the rows being the system states at t0, t1, t2, t3, t4); let's look at its time derivative:

In [23]:

gradient_t_at_x2 = np.gradient(f_of_t_at_x2, delta_t)
gradient_t_at_x2

Out[23]:

array([-9.0817816 , -9.03841995, -8.94751865, -8.84887524, -8.79777149])

The above is the rate of change of the concentration, as the diffusion proceeds, at the position x2¶

At time t2, the midpoint in the simulation, the value is:

In [24]:

gradient_t_at_x2[2]

Out[24]:

-8.947518648702157

C. All said and done, we have collected the time derivative and the 2nd spacial derivative, at the point (x2, t2).¶

Do those values satisfy the diffusion equation?? Does the time derivative indeed equal the (diffusion rate) x (the 2nd spacial derivative)? Let's see:

In [25]:

gradient_t_at_x2[2]

Out[25]:

-8.947518648702157

In [26]:

diffusion_rate * second_gradient_x_at_t2[2]

Out[26]:

-7.877903914825475

D. The 2 value indeed roughly match - considering the coarseness of the large spacial grid, and the coarseness of estimating the derivatives numerically.¶

We have evidence that the diffusion equation may indeed be satisfied by the values we obtained from our diffusion algorithm, in the proximity of the point (x2, t2)

E. Finally, instead of just scrutining the match at the point x2, let's do that for all the points in space at time t2,¶

WITH THE EXCEPTION of the outmost points (because the numeric estimation of the derivatives gets very crummy at the boundary).

Let's first re-visit all the data once again:

In [27]:

all_history

Out[27]:

array([[50.        , 55.87785252, 59.51056516, 59.51056516, 55.87785252,
        50.        , 44.12214748, 40.48943484, 40.48943484, 44.12214748],
       [50.14694631, 55.82172403, 59.41974735, 59.41974735, 55.82172403,
        50.        , 44.17827597, 40.58025265, 40.58025265, 44.03132966],
       [50.28881576, 55.76980517, 59.32979676, 59.32979676, 55.76613151,
        50.        , 44.23386849, 40.67020324, 40.66652958, 43.94505274],
       [50.42584049, 55.72178022, 59.24079697, 59.24070513, 55.71106985,
        50.        , 44.28893015, 40.75920303, 40.7485845 , 43.86308966],
       [50.55823898, 55.67735715, 59.15281926, 59.15246655, 55.65653399,
        50.        , 44.34346372, 40.84718074, 40.82671259, 43.78522703]])

In [28]:

# The following is just an expanded version of what we did before for the time derivative; 
# instead of just considering 1 column,
# like we did before, we're now repeating the computations along all columns
# (the computations are applied vertically, "along axis 0" in Numpy-speak)
gradient_t = np.apply_along_axis(np.gradient, 0, all_history, delta_t)
gradient_t

Out[28]:

array([[14.69463131, -5.61284971, -9.0817816 , -9.0817816 , -5.61284971,
         0.        ,  5.61284971,  9.0817816 ,  9.0817816 , -9.0817816 ],
       [14.44078779, -5.40236784, -9.03841995, -9.03841995, -5.58605073,
         0.        ,  5.58605073,  9.03841995,  8.85473706, -8.85473706],
       [13.9447089 , -4.99719025, -8.94751865, -8.95211072, -5.5327087 ,
         0.        ,  5.5327087 ,  8.94751865,  8.41659228, -8.41200021],
       [13.47116142, -4.62240099, -8.84887524, -8.86651087, -5.47987603,
         0.        ,  5.47976123,  8.84887524,  8.00915063, -7.99128539],
       [13.23984932, -4.44230744, -8.79777149, -8.8238586 , -5.45358643,
         0.        ,  5.45335682,  8.79777149,  7.81280921, -7.7862629 ]])

In [29]:

# Again, we focus on time t2 (the 3rd row), to stay away from the edges
gradient_t_at_t2 = gradient_t[2]
gradient_t_at_t2

Out[29]:

array([13.9447089 , -4.99719025, -8.94751865, -8.95211072, -5.5327087 ,
        0.        ,  5.5327087 ,  8.94751865,  8.41659228, -8.41200021])

Note the value -8.94751865, 3rd from left : that's the single value we looked at before (at point x2)

Time to again check the match of the two sides of the diffusion equation¶

In [30]:

lhs = gradient_t_at_t2   # The LEFT-hand side of the diffusion equation, as a vector for all the spacial points at time t2
lhs

Out[30]:

array([13.9447089 , -4.99719025, -8.94751865, -8.95211072, -5.5327087 ,
        0.        ,  5.5327087 ,  8.94751865,  8.41659228, -8.41200021])

In [31]:

rhs = diffusion_rate * second_gradient_x_at_t2  # The RIGHT-hand side of the diffusion equation, again as a vector
rhs

Out[31]:

array([-2.40124727, -4.62624201, -7.87790391, -8.05611773, -4.9803736 ,
        0.        ,  4.97807756,  7.87790391,  6.32774077,  4.10274602])

The left-hand side and the right-hand side of the diffusion equation appear to generally agree, except at the boundary points, where our approximations are just too crummy¶

In [32]:

lhs - rhs

Out[32]:

array([ 16.34595617,  -0.37094824,  -1.06961473,  -0.89599299,
        -0.5523351 ,   0.        ,   0.55463113,   1.06961473,
         2.08885151, -12.51474623])

In [33]:

# Here we use a handy function to compare two equal-sized vectors,
# while opting to disregarding a specified number of entries at each edge.
# It returns the Euclidean distance ("L2 norm") of the shortened vectors
Numerical.compare_vectors(lhs, rhs, trim_edges=1)

Out[33]:

2.8643581811669576

IMPORTANT: all values in this experiment are VERY coarse, because of the large effective delta_x (the tiny number of bins.)¶

In part2 (notebook "validate_diffusion_2"), much-better approximations will get looked at!

In [ ]: