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

268 8 SolvingOrdinaryDifferentialEquations whichdiffers fromu1 in (8.78)byanamount 12Δt 2ω2u0. Because of the equivalenceof (8.76)with the Euler-Cromerscheme, the numer- ical results will have the same nice properties such as a constant amplitude. There will be a phase error as in the Euler-Cromer scheme, but this error is effectively reducedbyreducingΔt, as alreadydemonstrated. The implementation of (8.78) and (8.76) is straightforward in a function (file osc_2nd_order.py): import numpy as np def osc_2nd_order(U_0, omega, dt, T): """ Solve u’’ + omega**2*u = 0 for t in (0,T], u(0)=U_0 and u’(0)=0, by a central finite difference method with time step dt. """ Nt = int(round(T/dt)) u = np.zeros(Nt+1) t = np.linspace(0, Nt*dt, Nt+1) u[0] = U_0 u[1] = u[0] - 0.5*dt**2*omega**2*u[0] for n in range(1, Nt): u[n+1] = 2*u[n] - u[n-1] - dt**2*omega**2*u[n] return u, t 8.4.13 AFiniteDifferenceMethod;LinearDamping A key issue is how to generalize the scheme from Sect. 8.4.12 to a differential equationwithmore terms. We start with the case of a lineardamping termf(u′)= bu′, a possiblynonlinear springforces(u), andanexcitationforceF(t): mu′′+bu′+s(u)=F(t), u(0)=U0, u′(0)=0, t∈ (0,T] . (8.79) We need to find the appropriate difference approximation tou′ in the bu′ term. A goodchoice is thecentereddifference u′(tn)≈ u n+1−un−1 2Δt . (8.80) Sampling theequationat a timepoint tn, mu′′(tn)+bu′(tn)+s(un)=F(tn), and inserting the finitedifferenceapproximations tou′′ andu′ results in m un+1−2un+un−1 Δt2 +bu n+1−un−1 2Δt +s(un)=Fn, (8.81)