Computational Physics Differential Equations

1. Computational PhysicsDifferential Equations Dr. Guy Tel-Zur Autumn Colors, by Bobby Mikul, http://www.publicdomainpictures.net Version 10-11-2010 18:30

2. Agenda • MHJ Chapter 13 & Koonin Chapter 2 • How to solve ODE using Matlab • Scilab

3. Topics • Defining the scope of the discussion • Simple methods • Multi-Step methods • Runge-Kutta • Case Studies - Pendulum

4. The scope of the discussion For a higher order ODE  a set of coupled 1st order equations:

5. Simple methods Euler method: Integration using higher order accuracy: Taylor series expansion:

6. Local error! Better than Euler’s method but useful only when it is easy to differentiate f(x,y)

7. An Example Boundary conditiion Let’s solve:

9. FORTRAN code - I Euler’s method Demo: virtual box, folder: ~/fortran, program: chap2a.f Compilation: telzur@linux1:~/fortran$ fort77 -o chap2a chap2a.f

10. FORTRAN code - II Taylor’s series method FUNC(X,Y)=-X*Y C-------scientific computing course, lecture 05 - differential equations C Guy Tel-Zur, April 2011 C example chap2a from Koonin C compile under ubuntu using: fort77 -o chap2a_taylor chap2a_taylor.f 20 PRINT *,'ENTER STEP SIZE' READ *,h IF (H.LE.0) STOP NSTEP=3./h Y=1. DO 10 IX=0,NSTEP-1 X=IX*H Y=Y+H*FUNC(X,Y)+0.5*H*H*(-Y-FUNC(X,Y)*X) DIFF=EXP(-0.5*(X+H)**2)-Y PRINT *,IX,X+H,Y,DIFF 10 CONTINUE GOTO 20 END

11. The Results Taylor’s Euler’s IX X Y DIFF Y DIFF 0 0.5 1. -.117503099 .875 .00749690272 1 1.0 .75 -.143469334 .57421875 .0323119089 2 1.5 .375 -.0503475331 .287109375 .0375430919 3 2.0 .09375 .0415852815 .116638184 .0186970998 4 2.5 0. .0439369343 .0437393188 .000197614776 5 3.0 0. .0111089963 .0177690983 -.00666010194

12. Multi-Step methodsAdams-Bashforth 2 steps: 4 steps:

13. (So far) Explicit methods Future = Function(Present && Past) Implicit methods Future = Function(Future && Present && Past)

14. Let’s calculate dy/dx at a mid way between lattice points: Let’s replace: Rearrange: This is a recursion relation!

15. A simplification occurs if f(x,y)=y*g(x), then the recursion equation becomes: An example, suppose g(x)=-x yn+1=(1-xnh/2)/(1+xn+1h/2)yn This can be easily calculated, for example: Calculate y(x=1) for h=0.5 X0=0, y(0)=1 X1=0.5, y(0.5)=? x2=1.0, y(1.0)=?

16. The solution: Error=-0.01569

17. Predictor-Corrector method

18. Predictor-Corrector method

19. Runge-Kutta

22. Proceed to: Physics examples: Ideal harmonic oscillator – section 13.6.1

23. Physics Project – The pendulum, 13.7 I use a modified the C++ code from: http://www.fys.uio.no/compphys/cp/programs/FYS3150/chapter13/cpp/program2.cpp fout.close fout.close() Demo on folder: \Lectures\05\CPP

24. ODEs in Matlab function dydt = odefun(t,y) a=0.001; b=1.0; dydt=b*t*sin(t)+a*t*t; Usage: [t1, y1]=ode23(@odefun,[0 100],0); [t2, y2]=ode45(@odefun,[0 100],0); plot(t1,y1,’r’); hold on plot(t2,y2,’b’); hold off Demo folder: C:\Users\telzur\Documents\Weizmann\ScientificComputing\SC2011B\Lectures\05\Matlab

25. Output

26. http://www.scilab.org

31. Parallel tools for Multi-Core and Distributed Parallel Computing

32. In preparation Parallel execution A new function (parallel_run) allows parallel computations and leverages multicore architectures and their capacities. In future Scilab versions: http://help.scilab.org/docs/5.3.1/en_US/parallel_run.html

38. Xcos demo