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

242 8 SolvingOrdinaryDifferentialEquations resulting in thecomputationalscheme un+1 =un+Δtvn, (8.47) vn+1 =vn−Δtω2un. (8.48) 8.4.3 ProgrammingtheFEScheme;theSpecialCase Asimpleprogramfor (8.47)–(8.48)follows thesame ideasas inSect. 8.3.3: import numpy as np import matplotlib.pyplot as plt omega = 2 P = 2*np.pi/omega dt = P/20 T = 3*P N_t = int(round(T/dt)) t = np.linspace(0, N_t*dt, N_t+1) u = np.zeros(N_t+1) v = np.zeros(N_t+1) # Initial condition X_0 = 2 u[0] = X_0 v[0] = 0 # Step equations forward in time for n in range(N_t): u[n+1] = u[n] + dt*v[n] v[n+1] = v[n] - dt*omega**2*u[n] fig = plt.figure() l1, l2 = plt.plot(t, u, ’b-’, t, X_0*np.cos(omega*t), ’r--’) fig.legend((l1, l2), (’numerical’, ’exact’), ’upper right’) plt.xlabel(’t’) plt.savefig(’tmp.pdf’); plt.savefig(’tmp.png’) plt.show() (Seefileosc_FE.py.) Sincewealreadyknowtheexactsolutionasu(t)=X0cosωt,wehavereasoned asfollowstofindanappropriatesimulationinterval[0,T]andalsohowmanypoints we should choose.The solution hasa periodP = 2π/ω. (The periodP is the time difference between two peaks of the u(t) ∼ cosωt curve.) Simulating for three periods of the cosine function,T = 3P, and choosingΔt such that there are 20 intervalsper periodgivesΔt =P/20 and a total ofNt =T/Δt intervals. The rest of the programisastraightforwardcodingof theForwardEuler scheme. Figure 8.19 shows a comparison between the numerical solution and the exact solution of the differential equation. To our surprise, the numerical solution looks wrong. Is this discrepancy due to a programming error or a problem with the ForwardEulermethod?