ODE placement in an ODE system give different results

Question

I have the following code:

def multirk4(funcs, x0, y0, step, xmax):
    n = len(funcs)
    table = [[x0] + y0]
    f1, f2, f3, f4 = [0]*n, [0]*n, [0]*n, [0]*n
    while x0 < xmax:
        y1 = [0]*n
        for i in range(n): f1[i] = funcs[i](x0, y0)
        for j in range(n): y1[j] = y0[j] + (0.5*step*f1[j])
        for i in range(n): f2[i] = funcs[i]((x0+(0.5*step)), y1)
        for j in range(n): y1[j] = y0[j] + (0.5*step*f2[j])
        for i in range(n): f3[i] = funcs[i]((x0+(0.5*step)), y1)
        for j in range(n): y1[j] = y0[j] + (step*f3[j])
        for i in range(n): f4[i] = funcs[i]((x0+step), y1)
        x0 = x0 + step
        for i in range(n): y1[i] = y0[i] + (step * \
            (f1[i] + (2.0*f2[i]) + (2.0*f3[i]) + f4[i]) / 6.0)
        table.append([x0] + y1)
        y0 = y1
    return table

system1 = range(2)
system2 = range(2)
y = range(2)

y[0] = 0.0
y[1] = 0.0

def mRNA(t, y): return 6e-8 - (0.01 * y[0])
def protein(t, y): return (0.05 * y[0]) - (0.001 * y[1])

system1[0] = mRNA
system1[1] = protein

system2[0] = protein
system2[1] = mRNA

t0 = 0.0
tmax = 100.0
dt = 0.1

s1 = multirk4(system1, t0, y, dt, tmax)
s2 = multirk4(system2, t0, y, dt, tmax)

for i in range(len(s1)):
    print ','.join([str(s1[i][0]), str(s1[i][1]), str(s1[i][2]), 
                    str(s2[i][1]), str(s2[i][2])])

I found that s1 and s2 give different results, which is essentially a different ordering of the equations (mRNA and protein).

Is there any amendments that I can make so that the ordering does not make a difference?

Thanks in advance.

Are you using Runge kutta 4? – kolonel Nov 01 '14 at 02:25 — kolonel, Nov 01 '14 at 02:25
Yes, I am using Runge Kutta 4 – Maurice Ling Nov 02 '14 at 08:39 — Maurice Ling, Nov 02 '14 at 08:39

Simon · Answer 1 · 2014-11-01T03:03:43.347

The problem seems to be that the functions are being swapped between indices 0 and 1 between cases s1 and s2 but the y values used by those functions are not being swapped. If we swap the y values as well, we get a consistent result. We can do this by replacing the function definitions for mRNA() and protein() and the setup code for the system1[] and system2[] arrays with the following code:

def mRNA1(t, y):
    print t,"\t",y[0]
    return 6e-8 - (0.01 * y[0])

def protein1(t, y): 
    return (0.05 * y[0]) - (0.001 * y[1])

def mRNA2(t, y):
    print t,"\t",y[0]
    return 6e-8 - (0.01 * y[1])

def protein2(t, y): 
    return (0.05 * y[1]) - (0.001 * y[0])

system1[0] = mRNA1
system1[1] = protein1

system2[0] = protein2
system2[1] = mRNA2

... we get output for which the last line is:

100.1,3.79492952635e-06,1.06680785809e-05,1.06680785809e-05,3.79492952635e-06

... as we would expect.

Parenthetically, I'd suggest using scipy.integrate.ode instead of writing custom code to solve the ODEs if the intention is to get a good solution to the ODEs rather than to learn how to write ODE solver code. All the solvers available in scipy will give a better result than using a 4th order Runge Kutta method.

Thanks, Simon. I am still puzzled about the difference. Can you kindly explain that to me or point me somewhere? Thanks a lot. — Maurice Ling, Nov 02 '14 at 08:40
The result of computing `mRNA()` is going into `system1[0]` (hence it was calculating the value of `y[0]`) in the first case and `system2[1]` (hence calculating the value of `y[1]`) in the second case. However in the original code it was always using the value of `y[0]` as input to compute the `y` value at the next timestep, regardless of whether its output was going into `y[0]` or `y[1]`. Conversely, the same principle applied to the `protein()` function. Swapping the `y[]` indices as well as the `system[]` indices makes them consistent and therefore produces the expected output. — Simon, Nov 02 '14 at 09:49

ODE placement in an ODE system give different results

1 Answers1