1
$\begingroup$

I strongly suspect that matrix methods are the right answer here, but I am having trouble visualizing how to set it up. Imagine a large number of coupled functions:

AB[m,n][t]

where m and n are whole number indices that can range up to hundreds or low thousands, and t is a continuous time variable. Think of AB[m,n] as the concentration of AB of type [m,n], and these concentrations of the different types can evolve with time. In addition, we have a single extra function:

B[t]

that is also a concentration that evolves with time, but for which there is only a single type. Initial conditions are:

AB[a,0][0] = A0
AB[anythingelse][0] = 0
B[0] = B0

Here, a is a constant whole number in the hundreds to low thousands.

AB[m,n] can transform into other types by two processes:

AB[m,n] + B --> AB[m-1,n+b]
AB[m,n] --> AB[m-1,n-1]

Here, b is a constant whole number significantly less than a. This means that the differential equations governing the evolution of the system are:

D[AB[m, n][t], t] == -k1[m, n] AB[m, n][t] B[t] 
                     +k1[m + 1, n - b] AB[m + 1, n - b][t] B[t] 
                     -k2[m, n] AB[m, n][t] 
                     +k2[m + 1, n + 1] AB[m + 1, n + 1][t]
D[B[t], t] == -Sum[k1[m, n] AB[m, n][t],{m,0,a},{n,0,Infinity}] B[t]

Good guesses for the forms for k1 and k2, which we can use for testing, are:

k1[m_] = (m/a) k3
k2[m_, n_] = k4 n/(n + k5/m) + n k6 k1[m]

Here, k3, k4, k5, and k6 are positive real numbers.

How on Earth do I organize this system in a way to numerically solve this system using NDSolve (or ParametricNDSolve) or similar? The ultimate goal will be to have a measurable function that would be something like:

Sum[m AB[m,n],{m,0,a},{n,0,infinity}]

This function would then be fit to experimental data, varying k3-k6, and maybe a and b. Later generalizations may be to have a broader distribution of starting concentrations rather than all being identical, a range of different possible bs, etc.

$\endgroup$
  • $\begingroup$ Are you really trying to solve a system of ODEs with an infinite number of variables (that is, does n really range from 0 to Infinity)? If not, have you tried creating a matrix equation to represent the ODEs, something like: AB'[t] == m1 . AB[t] . m2 B[t] + n1 . AB[t] . n2, where m1, m2, n1 and n2 encode the information about your k variables. $\endgroup$ – Carl Woll Nov 11 '19 at 22:43
  • $\begingroup$ Given the mechanism, the maximum value that n could take is a times b. I will look at your suggestion. $\endgroup$ – Kevin Ausman Nov 11 '19 at 22:55
  • $\begingroup$ Looks like I found a way to make your suggestion work! Thank you! I will post a completed answer shortly, once I debug one aspect. $\endgroup$ – Kevin Ausman Nov 12 '19 at 0:58
3
$\begingroup$

Ok, here we go. Let's start with a relatively simple system:

a = 20;
b = 4;

Also, my equations above are slightly off. The first transformation should have been:

AB[m,n] + B --> AB[m-1,n+(b-1)]

which changes the differential equations somewhat. First, let's put in our definitions for the rate constants.

k1[m_, n_] = (m/a) k3
k2[m_, n_] = k4 n m/(n m + k5) + n k6 k1[m]

And some initial guess values:

k3 = k4 = k5 = k6 = 1;
AB0 = 10;
B0 = 50;

Now we'll make a vector for all of the species:

ABvector = Flatten@Table[AB[m, n][t], {m, 0, a}, {n, 0, a (b - 1)}];

We can index AB[m,n] via:

index[m_, n_] = m (a (b - 1) + 1) + n + 1;
AB[[index[m,n]]]

Now the k1 and k2 matrices:

k1matrix = SparseArray[
             Flatten@Table[{{index[m, n], index[m, n]} -> -k1[m, n], 
                            {index[m - 1, n + (b - 1)], index[m, n]} -> k1[m, n]}, 
                           {m, 1, a}, 
                           {n, 0, a (b - 1) - (b - 1)}], 
             {(a + 1) (1 + a (b - 1)), (a + 1) (1 + a (b - 1))}];
k2matrix = SparseArray[
             Flatten@Table[{{index[m, n], index[m, n]} -> -k2[m, n], 
                            {index[m - 1, n - 1], index[m, n]} -> k1[m, n]}, 
                           {m, 1, a}, 
                           {n, 1, a (b - 1)}], 
             {(a + 1) (1 + a (b - 1)), (a + 1) (1 + a (b - 1))}];

We can now define all of our rate equations (except for B) like this:

rateEq = Thread[D[ABvector, t] == k1matrix.ABvector B[t] + k2matrix.ABvector];

And our initial conditions are (except for B):

initCond = Thread[(ABvector /. t -> 0) == SparseArray[{index[a, 0] -> AB0}, 
                                           {(a + 1) (1 + a (b - 1))}]];

We can now run NDSolve:

result = NDSolve[Join[rateEq, initCond, 
                      {B[0] == B0, B'[t] == -Sum[k1[m, n] B[t] AB[m, n][t], 
                                                 {m, 1, a}, 
                                                 {n, 0, a (b - 1) - (b - 1)}]}], 
                 Append[ABvector, B[t]], {t, 0, 10}, WorkingPrecision -> 30];

And we can visualize the results like so:

fn[t_] = Table[AB[m, n][t], {m, 0, a}, {n, 0, a (b - 1)}] /. (result[[1]]);
frameTable = Table[ListPointPlot3D[fn[t], ImageSize -> Large, Filling -> 0, 
                   FillingStyle -> Thickness[0.01], ColorFunction -> "Rainbow", 
                   AxesLabel -> {"n", "m", "Concentration"}, 
                   ViewPoint -> {1.3, -2.4, 1.}], {t, 0, 3, 0.1}];
ListAnimate[frameTable]

enter image description here

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.