Page - 208 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python 3.6, Volume Second Edition

208 8 SolvingOrdinaryDifferentialEquations errors are “small”, also the second straight line segment shouldbe close to the true solutioncurve. 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 errormay grow with thenumberof time steps, but since the total time interval is not too large, andsince wemaychoosea very small timestep onmoderncomputers, this couldstill 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.1s for a first try. A simple versionof 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, somemoreflexibility is introducedhere,usingrangeandN in thefor loopheader. Runningthecodegives theplot shownin Fig.8.2. Thislookspromising!Notsurprisingly,theerrorgrowswithtime,reachingabout 2.64L 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.30L, while an N value of 3 · 106 gives an error of 3 · 10−5L. 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 changebo- to, e.g., onlyb.