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

4.3 OscillatingOne-DimensionalSystems 127 # 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*cos(omega*t), ’r--’) fig.legend((l1, l2), (’numerical’, ’exact’), ’upper left’) plt.xlabel(’t’) plt.show() plt.savefig(’tmp.pdf’); plt.savefig(’tmp.png’) (Seefileosc_FE.py.) Sincewealreadyknowtheexactsolutionasu.t/DX0cos!t,wehavereasoned asfollowstofindanappropriatesimulationinterval Œ0;T andalsohowmanypoints weshouldchoose. ThesolutionhasaperiodP D2 =!. (TheperiodP is the time difference between two peaks of theu.t/ cos!t curve.) Simulating for three periods of the cosine function,T D 3P , and choosing t such that there are 20 intervalsperperiodgives t DP=20anda totalofNt DT= t intervals. Therest of theprogramis a straightforwardcodingof theForwardEuler scheme. Figure 4.16 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 discrepancydue to a programmingerror or a problemwith theFor- wardEulermethod? First of all, even before trying to run the program, you should sit down and compute twosteps in the timeloopwithacalculator soyouhavesomeintermediate results to comparewith. UsingX0 D 2,dt D 0:157079632679, and! D 2, we Fig.4.16 Simulationof anoscillatingsystem