Runge-Kutta 4 plot time range not limited by tmax

Question

I am trying to solve an ODE using the 4th order Runge Kutta method. The problem is that the defined tmax in my code does not seem to affect the time range of the plot. Changing tmax only affects the curve of the plot. Instead the time range is dependent on nT.

#generated electricity
g = 1.0 

#rate constants
k_diss = 0.1
k_em = .01

#4th order RK
def RK4(func, yi, ti, dt):
    k1 = dt*func(yi, ti)
    k2 = dt*func(yi + 0.5 * k1, ti + 0.5 * dt)
    k3 = dt*func(yi + 0.5 * k2, ti + 0.5 * dt)
    k4 = dt*func(yi + k3,ti + dt)
    return yi + 1./6*(k1 + 2*k2 + 2*k3 + k4)

#variation of cEX
def df(t, c, k_diss, k_em):
    return np.array([g - 2*k_diss*c[0] - k_em*c[0]])

#time
tmax = 30.0
nT = 500
T = np.linspace(0,tmax,nT)
dT = tmax/(nT-1)

#initial condition
C0, T0 = 1.0, 0.0

#storage
C = np.zeros(nT)
C[0] = C0
 
#ODE solver

r = ode(df).set_integrator('dopri5')
r.set_initial_value(C0,T0).set_f_params(k_diss,k_em)

#loop over time
iT = 1

for i in range(1,nT):
    C[iT] = r.integrate(r.t + dT)
    iT += 1

On top of that the plot should be decreasing (can be seen below), not increasing. I have based my code on the following (there tmax does operate the way it is supposed to):

def RK4(func, yi, ti, dt):
    k1 = dt*func(yi,ti)
    k2 = dt*func(yi+0.5*k1,ti+0.5*dt)
    k3 = dt*func(yi+0.5*k2,ti+0.5*dt)
    k4 = dt*func(yi+k3,ti+dt)
    return yi + 1./6*  ( k1+2*k2+2*k3+k4 )

# funtion that defines the concentration of CA
def df(c,t):
    return np.array([-k*c[0],k*c[0]])   # return the value

# time argument
tmax,nT = 30.0,50
T = np.linspace(0,tmax,nT)
dT = tmax/(nT-1)

# rate constant
k = 0.2

# declare the array to store
RK = np.zeros((2,nT))

# initial condition
RK[0,0] = 1.

for iT in range(1,nT):
    RK[:,iT] = RK4(df,RK[:,iT-1],T[iT-1],dT)

The plot:

Honestly, your first code is a mess. You define a lot of things that you do not use (including `T`). Also, what do you mean by “time range of the plot”? There is no plot in your code. You get a one-dimensional array. Its length is obviously `nT` (because that’s how you initialise it). Please read [mcve] and [edit] your question accordingly. — Wrzlprmft, Apr 08 '17 at 15:11
@Peaceful: I fail to make sense of your request. `func` is internal variable of the `RK4` function (which in turn is irrelevant, as it is never used again). — Wrzlprmft, Apr 08 '17 at 15:15

score 0 · Answer 1 · answered Apr 27 '17 at 08:05

I changed the end of the first script to

#loop over time
for iT in range(1,nT):
    C[iT] = r.integrate(T[iT])

import matplotlib.pyplot as plt
plt.plot(T,C)
plt.xlabel("time"); plt.ylabel("cEX");
plt.show()

The first change since it is actually your intent that C[i] is the value at time T[i], and as there is a way to express this consistently, why not use it.

Since you did not add your plot commands, one can only guess your error. In all probability it was that you did not include the time in your plot command. If you only write plt.plot(C), then the horizontal axis is filled with the integer indices of C, since the plot procedure is not given knowledge about any other time scale.

Runge-Kutta 4 plot time range not limited by tmax

1 Answers1