Computational Physics using MATLAB®

1 kevin Berwick Page 1 Computational Physics using MATLAB kevin Berwick Page 2 Table of Contents 6 1. Uranium Decay .. 7 3. The Pendulum .. 9 Solution using the Euler method .. 9 Solution using the Euler-Cromer method.. 10 Simple Harmonic motion example using a variety of numerical approaches .. 11 Solution for a damped pendulum using the Euler-Cromer method.. 16 Solution for a non-linear, damped, driven pendulum :- the Physical pendulum, using the Euler-Cromer method.. 18 Bifurcation diagram for the pendulum .. 24 The Lorenz Model .. 26 4. The Solar System .. 28 Kepler s Laws .. 28 Ex Planetary motion results using different time steps ..30 Orbits using different force laws .. 35 Precession of the perihelion of Mercury.. 40 The three body problem and the effect of Jupiter on Earth .. 48 Chaotic tumbling of Hyperion .. 53 5. Potentials and Fields .. 60 Solution of Laplace s equation using the Jacobi relaxation method.

2 60 Solution of Laplace s equation for a hollow metallic prism with a solid, metallic inner conductor.. 63 Solution of Laplace s equation for a finite sized capacitor .. 66 Exercise and the Successive Over Relaxation Algorithm .. 70 Potentials and fields near Electric charges, Poisson s Equation .. 75 6. Waves .. 78 Waves on a string .. 78 Waves on a string with free ends .. 81 Frequency spectrum of waves on a string .. 83 7. Random Systems .. 87 Random walk simulation .. 87 Random walk simulation with random path lengths..89 10. Quantum Mechanics .. 91 Time independent Schrodinger equation. Shooting method.. 91 kevin Berwick Page 3 Wavepacket construction .. 93 Time Dependent Schrodinger equation in One dimension. Leapfrog method.. 95 Time Dependent Schrodinger equation in two dimensions. Leapfrog method.. 99 kevin Berwick Page 4 Table of Figures Figure 1. Uranium decay as a function of time.

3 8 Figure 2. Simple Pendulum - Euler Method .. 9 Figure 3. Simple Pendulum: Euler - Cromer method .. 10 Figure 4. Simple pendulum solution using Euler, Euler Cromer, Runge Kutta and Matlab ODE45 solver.. 15 Figure 5. The damped pendulum using the Euler-Cromer method .. 17 Figure 6. Results from Physical pendulum, using the Euler-Cromer method, F_drive = 19 Figure from Physical pendulum, using the Euler-Cromer method, F_drive = ..20 Figure 8. Results from Physical pendulum, using the Euler-Cromer method, F_drive = 21 Figure 9. Results from Physical pendulum, using the Euler-Cromer method, F_Drive= . 21 Figure 10. Increase resolution with npoints= from Physical pendulum, using the Euler-Cromer method, F_Drive= .. 22 Figure 11. Poincare section (Strange attractor) Omega as a function of theta. F_Drive = . 23 Figure 12. Bifurcation diagram for the pendulum .. 25 Figure 13. Variation of z as a function of time and corresponding strange attractor.

4 27 Figure 14. Simulation of Earth orbit around the Sun .. 29 Figure 15. Simulation of Earth orbit with time step of .. 31 Figure 16. Simulation of Earth orbit, initial y velocity of 4, time step is .. 32 Figure of Earth orbit, initial y velocity of 4, time step is .. 33 Figure 18. Simulation of Earth orbit, initial y velocity of 8, time step is 2500 points and Runge Kutta method .. 33 Figure for an initial y velocity of 8, dt is , npoints=2500. The Runge Kutta Method is used here .. 35 Figure 20. Orbit for a force law with =2. The time step is years.. 37 Figure 21. Orbit for a force law with = The time step is years.. 39 Figure 22. Orbit for a force law with =3.. 40 Figure 23. Orbit orientation as a function of time .. 45 Figure 24. Calculated precession rate of Mercury .. 47 Figure 25. Simulation of solar system containing Jupiter and Earth .. 50 Figure 26. Simulation of solar system containing Jupiter and Earth.

5 Jupiter mass is 10 X actual value.. 51 Figure of solar system containing Jupiter and Earth. Jupiter mass is 1000 X actual value, ignoring perturbation of the Sun.. 52 Figure of Hyperion. The initial velocity in the y direction was 1 HU/Hyperion year. This gave a circular orbit. Note from the results that the tumbling is not chaotic under these conditions.. 56 Figure of Hyperion. The initial velocity in the y direction was 5 HU/Hyperion year. This gave a circular orbit. Note from the results that the tumbling is chaotic under these conditions.. 59 Figure 30. Equipotential surface for geometry depicted in Figure in the book .. 62 Figure surface for hollow metallic prism with a solid metallic inner conductor held at V=1.. 66 kevin Berwick Page 5 Figure surface for a finite sized capacitor.. 69 Figure 33. Equipotential contours near a finite sized capacitor.. 69 Figure surface in region of a simple capacitor as calculated using the SOR code for a 60 X 60 grid.

6 The convergence criterion was that the simulation was halted when the difference in successively calculated surfaces was less than 10-5 per 73 Figure of iterations required for Jacobi method vs L for a simple capacitor. The convergence criterion was that the simulation was halted when the difference in successively calculated surfaces was less than 10-5 per site.. 74 Figure of iterations required for SOR method vs L for a simple capacitor. The convergence criterion was that the simulation was halted when the difference in successively calculated surfaces was less than 10-5 per site.. 74 Figure 37. Equipotential surface near a point charge at the center of a 20X20 metal box. The Jacobi relaxation method was used . The plot on the right compares the numerical and analytical (as obtained from Coulomb s Law).. 77 Figure 38. Waves propagating on a string with fixed ends .. 79 Figure 39. Signal from a vibrating string and Power spectrum.

7 Signal excited with Gaussian pluck centred at the middle of the string and the displacement 5% from the end of the string was recorded..86 Figure 40. x^2 as a function of step number. Step length = 1. Average of 500 walks. Also shown is a linear fit to the data.. 88 Figure 41. x^2 as a function of step number. Step length = random value betwen +/-1. Average of 500 walks. Also shown is a linear fit to the data.. 90 Figure 42. Calculated wavefunction using the shooting method. The wall(s) of the box are at x=(-)1. The value of Vo used was 1000 giving ground-state energy of Analytical value is Wavefunctions are not normalised.. 93 Figure 43. Composition of wavepacket. ko = 500, x0= , sigma^2= .. 94 Figure 44. Wavepacket reflection from potential cliff at x= The potential was V=0 for x< and V=-1e6 for x> Values used for initial wavepacket were x_0= ,C=10, sigma_squared=1e-3, k_0=500. Simulation used delta_x=1e-3, delta_t=5e-8.

8 Time progresses left to right..98 Figure 45. Wavepacket reflection from potential wall at x= The potential was V=0 for x< and V=1e6 for x> Values used for initial wavepacket were x_0= ,C=10, sigma_squared=1e-3, k_0=500. Simulation used delta_x=1e-3, delta_t=5e-8. Time progresses left to right..98 Figure reflection from potential cliff at x= The potential was V=0 for x< and V=-1e3 for x> Values used for initial wavepacket were x_0= , y_0= ,C=10, sigma_squared= , k_0=40. Simulation used delta_x= , delta_t= .. 102 Figure 47. Wavepacket reflection from potential wall at x= The potential was V=0 for x< and V=1e3 for x> Values used for initial wavepacket were x_0= , y_0= ,C=10, sigma_squared= , k_0=40. Simulation used delta_x= , delta_t= .. 103 kevin Berwick Page 6 Preface I came across the book, Computational Physics , in the library here in the Dublin Institute of Technology in early 2012. Although I was only looking for one, quite specific piece of information, I had a quick look at the Contents page and decided it was worth a more detailed examination.

9 I hadn t looked at using numerical methods since leaving College almost a quarter century ago. I cannot remember much attention being paid to the fact that this stuff was meant to be done on a computer, presumably since desktop computers were still a bit of a novelty back then. And while all the usual methods, Euler, Runge-Kutta and others were covered, we didn t cover applications in much depth at all. It is very difficult to anticipate what will trigger an individual s intellectual curiosity but this book certainly gripped me. The applications were particularly well chosen and interesting. Since then, I have been working through the exercises intermittently for my own interest and have documented my efforts in this book, still a work in progress. Coincidentally, I had started to use MATLAB for teaching several other subjects around this time. MATLAB allows you to develop mathematical models quickly, using powerful language constructs, and is used in almost every Engineering School on Earth.

10 MATLAB has a particular strength in data visualisation, making it ideal for use for implementing the algorithms in this book. The Dublin Institute of Technology has existing links with Purdue University since, together with UPC Barcelona, it offers a joint Master's Degree with Purdue in Sustainability, Technology and Innovation via the Atlantis Programme. I travelled to Purdue for two weeks in Autumn 2012 to accelerate the completion of this personal project. I would like to thank a number of people who assisted in the production of this book. The authors of Computational Physics , Nick Giordano and Hisao Nakanishi from the Department of Physics at Purdue must be first on the list. I would like to thank both of them sincerely for their interest, hospitality and many useful discussions while I was at Purdue. They provided lot of useful advice on the Physics , and their enthusiasm for the project when initially proposed was very encouraging.