208
8 Solving Ordinary Differential Equations
errors are “small”, also the second straight line segment should be close to the true
solution curve.
What About the Errors? We realize that, in this way, we can work our way all
along the total time interval. Immediately, we suspect that the error may grow with
the number of time steps, but since the total time interval is not too large, and since
we may choose a very small time step on modern computers, this could still work!
Implementation and Performance Let us write down the code, which by choice
gets very similar to the code in Case 1, and see how it performs. We realize that,
for our strategy to work, the time steps should not be too large. However, during
these initial investigations of ours, our aim is first and foremost to check out the
computational idea. So, we pick a time step Δt = 0.1 s for a first try. A simple
version of the code (rate_exponential.py) may then read:
import numpy as np
import matplotlib.pyplot as plt
a = 0.0; b = 3.0
# time interval
N = 30
# number of time steps
dt = (b - a)/N
# time step (s)
V = np.zeros(N+1)
# numerically computed volume (L)
V[0] = 1
# initial volume
for i in range(0, N, 1):
V[i+1] = V[i] + dt*V[i]
# ...r is V now
time_exact = np.linspace(a, b, 1000)
V_exact = np.exp(time_exact)
# make exact solution (for plotting)
time = np.linspace(0, 3, N+1)
plt.plot(time, V, ’bo-’, time_exact, V_exact, ’r’)
plt.title(’Case 2’)
plt.legend([’numerical’,’exact’], loc=’upper left’)
plt.xlabel(’t (s)’)
plt.ylabel(’V (L)’)
plt.show()
To plot the exact solution, we just picked 1000 points in time, which we consider
“large enough” to get a representative curve. Compared to the code for Case 1,
some more flexibility is introduced here, using range and N in the for loop header.
Running the code gives the plot shown in Fig. 8.2.
This looks promising! Not surprisingly, the error grows with time, reaching about
2.64 L at the end. However, the time step is not particularly small, so we should
expect much more accurate computations if Δt is reduced. We skip showing the
plots, 3 but if we increase N from 30 to 300, the maximum error drops to 0.30 L,
while an N value of 3 · 10 6 gives an error of 3 · 10 −5 L. It seems we are on to
something!
3 With smaller time steps, it becomes inappropriate to use filled circles on the graph for the
numerical values. Thus, in the plot command, one should change bo- to, e.g., only b.
8 Solving Ordinary Differential Equations
errors are “small”, also the second straight line segment should be close to the true
solution curve.
What About the Errors? We realize that, in this way, we can work our way all
along the total time interval. Immediately, we suspect that the error may grow with
the number of time steps, but since the total time interval is not too large, and since
we may choose a very small time step on modern computers, this could still work!
Implementation and Performance Let us write down the code, which by choice
gets very similar to the code in Case 1, and see how it performs. We realize that,
for our strategy to work, the time steps should not be too large. However, during
these initial investigations of ours, our aim is first and foremost to check out the
computational idea. So, we pick a time step Δt = 0.1 s for a first try. A simple
version of the code (rate_exponential.py) may then read:
import numpy as np
import matplotlib.pyplot as plt
a = 0.0; b = 3.0
# time interval
N = 30
# number of time steps
dt = (b - a)/N
# time step (s)
V = np.zeros(N+1)
# numerically computed volume (L)
V[0] = 1
# initial volume
for i in range(0, N, 1):
V[i+1] = V[i] + dt*V[i]
# ...r is V now
time_exact = np.linspace(a, b, 1000)
V_exact = np.exp(time_exact)
# make exact solution (for plotting)
time = np.linspace(0, 3, N+1)
plt.plot(time, V, ’bo-’, time_exact, V_exact, ’r’)
plt.title(’Case 2’)
plt.legend([’numerical’,’exact’], loc=’upper left’)
plt.xlabel(’t (s)’)
plt.ylabel(’V (L)’)
plt.show()
To plot the exact solution, we just picked 1000 points in time, which we consider
“large enough” to get a representative curve. Compared to the code for Case 1,
some more flexibility is introduced here, using range and N in the for loop header.
Running the code gives the plot shown in Fig. 8.2.
This looks promising! Not surprisingly, the error grows with time, reaching about
2.64 L at the end. However, the time step is not particularly small, so we should
expect much more accurate computations if Δt is reduced. We skip showing the
plots, 3 but if we increase N from 30 to 300, the maximum error drops to 0.30 L,
while an N value of 3 · 10 6 gives an error of 3 · 10 −5 L. It seems we are on to
something!
3 With smaller time steps, it becomes inappropriate to use filled circles on the graph for the
numerical values. Thus, in the plot command, one should change bo- to, e.g., only b.