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

260 8 SolvingOrdinaryDifferentialEquations Thefileosc_EC_general.pyhasafunctionEulerCromerthat implementsthis method: def EulerCromer(f, s, F, m, T, U_0, V_0, dt): import numpy as np N_t = int(round(T/dt)) print(’N_t:’, N_t) t = np.linspace(0, N_t*dt, N_t+1) u = np.zeros(N_t+1) v = np.zeros(N_t+1) # Initial condition u[0] = U_0 v[0] = V_0 # Step equations forward in time for n in range(N_t): v[n+1] = v[n] + dt*(1./m)*(F(t[n]) - f(v[n]) - s(u[n])) u[n+1] = u[n] + dt*v[n+1] return u, v, t The Fourth Order Runge-Kutta Method The RK4 method just evaluates the right-handside of theODEsystem, ( 1 m (F(t)−s(u)−f(v)),v) for known values ofu, v, and t, so the method is very simple to use regardless of howthe functionss(u)andf(v)arechosen. 8.4.9 IllustrationofLinearDamping We consider an engineering system with a linear spring, s(u)= kx, and a viscous damper, where the damping force is proportional to u′, f(u′) = bu′, for some constant b > 0. This choice may model the vertical spring system in a car (but engineersoften like to illustrate such a system by a horizontallymovingmass, like the one depicted in Fig. 8.28). We may choose simple values for the constants to illustrate basic effects of damping(and later excitations). Choosing the oscillations to be the simpleu(t) = cost function in the undamped case, we may setm = 1, k = 1, b = 0.3, U0 = 1, V0 = 0. The following function implements this case: def linear_damping(): import numpy as np b = 0.3 f = lambda v: b*v s = lambda u: k*u F = lambda t: 0 m = 1 k = 1