import numpy as np import matplotlib.pyplot as plt from gekko import GEKKO from scipy.integrate import odeint # number of time points n = 15 # final time tf = 7.0 # initial concentration Ca0 = 5.0 # constants k = 1 # 1/s # method #1 Analytical solution # Ca(t) = Ca(0) * exp(-k*t) t = np.linspace(0,tf,n) Ca_m1 = Ca0 * np.exp(-k*t) # method #2 Euler's method t2 = np.arange(0,tf,0.5) n = len(t2) Ca_m2 = np.empty_like(t2) Ca_m2[0] = Ca0 # kmol/m^3 for i in range(1,n): dt = t2[i] - t2[i-1] Ca_m2[i] = Ca_m2[i-1] - k * Ca_m2[i-1] * dt t3 = np.arange(0,tf,1.5) n = len(t3) Ca_m3 = np.empty_like(t3) Ca_m3[0] = Ca0 # kmol/m^3 for i in range(1,n): dt = t3[i] - t3[i-1] Ca_m3[i] = Ca_m3[i-1] - k * Ca_m3[i-1] * dt t4 = np.arange(0,tf,2.1) n = len(t4) Ca_m4 = np.empty_like(t4) Ca_m4[0] = Ca0 # kmol/m^3 for i in range(1,n): dt = t4[i] - t4[i-1] Ca_m4[i] = Ca_m4[i-1] - k * Ca_m4[i-1] * dt # plot results plt.figure(1) plt.plot(t,Ca_m1,'r-',label='Ca (Analytical)') plt.plot(t2,Ca_m2,'k.-',label='Ca (Euler 0.5)') plt.plot(t3,Ca_m3,'bo-',label='Ca (Euler 1.5)') plt.plot(t4,Ca_m4,'y--',label='Ca (Euler 2.1)') plt.xlabel('Time (sec)') plt.ylabel(r'\$C_a (kmol/m^3)\$') plt.legend(loc='best') plt.show()