Solving Ordinary Differential Equations: Example Problems
Try to solve these problems before watching the solutions in the screencasts.
Example Problem 1
Solve the two ordinary differential equations simultaneously using Polymath. dCA/dt = -k1CA dCB/dt = k1CA -k2CB2
Initial conditions: at t = 0, CA0 = 10, CB0 = 2
Polymath code
# Solve ODEs for reaction A –> B –> C in batch reactor
d(CA) / d(t) = -k1*CA #mass balance for A
CA(0) = 10
d(CB) / d(t) = k1*CA – k2*CB^2 #mass balance for B
CB(0) = 2
k1 = 0.45
k2 = 0.05
t(0) = 0
t(f) = 8
See the corresponding Polymath file.
Example Problem 2
Use Mathematica to solve the two differential equations simultaneously.
dCA/dt = -k1CACB
dCB/dt = -2k1CACB
Initial conditions: at t = 0, CA0 = 10, CB0 = 12, k1 = 0.1
Mathematica code
(*Batch reactor A + 2B –> C ; r = kcAcB *)
soln = NDSolve[{cA'[t] == -0.1*cA[t]*cB[t],
cB'[t] == -0.1*2*cA[t]*cB[t],
cA[0] == 10., cB[0] == 12.},
{cA, cB}, {t, 0., 5.}];
plot = Plot[{cA[t], cB[t]} /. soln, {t, 0, 5}]
See the corresponding Mathematica file.
Solve ODEs with MATLAB
This MATLAB code solves the same problem as solved above with Polymath.
function ODE example
clear, clc, format short g, format compact
tspan = [0 8.]; % Range for the independent variable
y0 = [10.; 2.]; % Initial values for the dependent variables
disp(‘ Variable values at the initial point ‘);
disp([‘ t = ‘ num2str(tspan(1))]);
disp(‘ y dy/dt ‘);
disp([y0 ODEfun(tspan(1),y0)]);
[t,y]=ode45(@ODEfun,tspan,y0);
for i=1:size(y,2)
disp([‘ Solution for dependent variable y’ int2str(i)]);
disp([‘ t y’ int2str(i)]);
disp([t y(:,i)]);
plot(t,y(:,i));
title([‘ Plot of dependent variable y’ int2str(i)]);
xlabel(‘ Independent variable (t)’);
ylabel([‘ Dependent variable y’ int2str(i)]);
pause
end
function dYfuncvecdt = ODEfun(t,Yfuncvec);
CA = Yfuncvec(1);
CB = Yfuncvec(2);
k1 = 0.45;
k2 = 0.05;
%mass balance for A
dCAdt = 0 – (k1 * CA);
%mass balance for B
dCBdt = k1 * CA – (k2 * CB ^ 2);
dYfuncvecdt = [dCAdt; dCBdt];
Solve ODEs with Python
This installation guide may be useful if you are new Python user:
http://www.umich.edu/~elements/5e/tutorials/python_installation_tutorial.pdf
See the Python example file.
Python code
# The following packages must exist in your Python environment in order
# to plot solutions to differential equations in Python:
# – scipy
# – matplotlib
# The following 27 lines of code attempt to open these packages, and download them if you do not already have them.
# Note that if you use Anaconda or a similar Python distribution, then you likely already have these packages and
# would replace all of this code with these four lines:
# from scipy.integrate import odeint
# import numpy as np
# import matplotlib.pyplot as plt
# import matplotlib
def install(package):
subprocess.call([sys.executable, “-m”, “pip”, “install”, package])
packagesInstalled = False
if sys.version_info[0] < 3:
raise Exception(“\nYou are using a deprecated version of Python. Visit the Python website to update to Python 3.\n”)
try:
import numpy as np
from scipy.integrate import odeint
except ImportError:
print(“\nYou are missing Python packages necessary to run this file. This installation will only occur once.\n”)
packagesInstalled = True
install(“scipy”)
import numpy as np
from scipy.integrate import odeint
try:
import matplotlib.pyplot as plt
import matplotlib
except ImportError:
if not packagesInstalled:
print(“\nYou are missing Python packages necessary to run this file. This installation will only occur once.\n”)
install(“matplotlib”)
import matplotlib.pyplot as plt
import matplotlib
# —– Start of differential equation solver —– #
# declare constants
k1 = 0.45
k2 = 0.05
# initial conditions
cA0 = 10
cB0 = 0
# returns a 2D array, [cA'(t), cB'(t)]. This needs to be a function because
# the first argument of the proceeding ODE integrator must be a callable function.
def slope(y, t, k1, k2):
cA, cB = y
dydt = [-k1 * cA, k1 * cA – k2 * cB ** 2]
return dydt
# initial time
t0 = 0
# final time
tf = 60
# 1000 numbers equally spaced between t0 and tf
t = np.linspace(t0, tf, 1000)
# initial conditions
y0 = [cA0, cB0]
# odeint is an ordinary differential equation solver, part of the scipy library. See documentation at
# docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html
sol = odeint(slope, y0, t, args=(k1, k2))
# update font size on the plot
matplotlib.rcParams.update({‘font.size’: 22, ‘font.weight’: ‘normal’})
# this plot will be 12 inches wide, 8 inches high. Actual size may vary.
plt.figure(figsize=(12, 8))
# adds cA and cB to the plot in blue (‘b’) and green (‘g’)
plt.plot(t, sol[:, 0], ‘b’, label=’cA(t)’)
plt.plot(t, sol[:, 1], ‘g’, label=’cB(t)’)
# axis labels
plt.xlabel(‘t’)
plt.ylabel(‘concentration (mol/L)’)
# adds a grid, and finally, displays the plot in a new window.
plt.grid()
plt.show()