Page - 151 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

4.3 OscillatingOne-DimensionalSystems 151 first timelevel(4.78)implementstheinitialconditionu0.0/slightlymoreaccurately thanwhat is naturallydone in theEuler-Cromerscheme. The latterwill do v1 Dv0 t!2u0; u1 Du0C tv1 Du0 t2!2u0; whichdiffers fromu1 in (4.78)byanamount 1 2 t2!2u0. Becauseof theequivalenceof (4.76)with theEuler-Cromerscheme, thenumer- ical resultswill have the sameniceproperties 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 (4.78) and (4.76) is straightforward in a function (file osc_2nd_order.py): from numpy import zeros, linspace 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. """ dt = float(dt) Nt = int(round(T/dt)) u = zeros(Nt+1) t = 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 4.3.13 AFiniteDifferenceMethod;LinearDamping A key issue is how to generalize the scheme from Sect. 4.3.12 to a differential equationwithmore terms.Westartwith thecaseofa lineardampingtermf.u0/D bu0, apossiblynonlinearspring forces.u/, andanexcitation forceF.t/: mu00Cbu0Cs.u/DF.t/; u.0/DU0; u0.0/D0; t 2 .0;T : (4.79) We need to find the appropriate difference approximation to u0 in the bu0 term. Agoodchoice is thecentered difference u0.tn/ u nC1 un 1 2 t : (4.80) Sampling theequationat a timepoint tn, mu00.tn/Cbu0.tn/Cs.un/DF.tn/;