Download Implementing 4th Order Non-Adaptive & Adaptive Runge-Kutta Methods for ODEs - Prof. Ludger and more Assignments Physics Fundamentals in PDF only on Docsity!
I. DUE DATE : April 13^ Homework 7 , 10:3 0 A^ - M^ Methods in Theoretical Physics (Physics 5350) (PLEASE SHOW YOUR WORK) Runge 1) Non Using the outlines from class, write a 4 - -Kutta and Adaptive Runge KuttaAdaptive 4th (^) – order Runge- th Kutta (^) order Runge-Kutta program. This program should not yet be adaptive and should consist of a main program, the Runge subroutine to calculate the solution after each time step. derivatives as given by the ODE. The main program should record the-Kutta subroutine ( the solver), and a The user should only have access to the main program and the derivative subroutine. So think about what needs to be passed from one routine to the next. This program will be used in Part 3 below (so you might want to look at that problem now). 2) Write a program that uses the adaptive 4 Organize your program into a Adaptive Runge-Kutta Method main program and ath (^) order Runge stepper, solver, - Kutta method as discussed in class. and derivative subroutine. The user should only have access to the The stepper subroutine should get as inputs the initial values, x, an array of dependent variables, y, main program and the derivative subroutine. an initial stepsize h, and the required accuracy (use fractional accuracy). The stepper subroutine should retu variables ynew, and the adaptively adjusted stepsize hrn the updated values, xnew. new, an updated array of dependent The main program should record the solution after each call to the stepper subroutine. 3) Test your program with a problem for which you know the solution. It is up to you to either come up with an example of your own or use the orbit of a comet (as in Homework 6 and again outlined below) or use an even simpler example like the one from Lecture #27 (y=sin(sqrt(k) * t**2) ). The idea here is to test your program before using it for a problem for which you don’t know the answer (like problem 4). 4) Use the programs to study a chaotic system (the Lorenz model) (See attachment). Hand in your work and email your source code of the programs, f [email protected] unctions, and subroutines to
Example Outline for Problem 3) Kepler’s Laws (very similar to homework #6) Let’s consider just one planet in our solar system. According to Newton force is given by ’s law of gravitation this ! where Ms and Me are the masses of the Sun and the planet, r the distance between them, and G is^ Fg^ =^ GM^ rs^2^ Me^ , the gravitational constant. Let’s further assume that the position of the sun is fixed. calculate the position of the planet (or comet) as a function of time. Our goal is to If we express Newton’s 2nd^ law in Cartesian coordinates we have
!
d dt^2 x 2 = F Mg , ex !
d dt^2 y 2 = F Mg , ey where F From Figure 1 we findg,x and Fg,y are the x and y components of the gravitational force.
! with a similar result for F sun, which is located at the origin of the coordinate system.^ Fg g,y, x^. Here the negative sign reminds us that the force is directed toward the=^ "^ GM^ rs^2^ M^ e cos^ #^ =^ "^ GM^ rs^3^ M^ ex Following our usual we get approach and write each of the second-order ODE as two first-order ODEs
!
dv dtx = " GM r 3 s x !
dx dt = vx !
dv dty = " GM r 3 s y !
dy dt = vy
As a first step you need to specify the initial conditions. Next you need to find an appropriate time step, different values of Δt. Keeping in mind that on the year Δt. You can do this by runnin-mark the orbit of the Earth should return tog your program with its initial condition, you can determine an ‘acceptable’ value for your time step. Document your results. Part b) Elliptical Orbits So far, we have used a circular orbit. This is a good choice to test your program. However the more interesting orbits are those that are elliptical. For these orbits, we will use our adaptive Runge-Kutta scheme. Initially, test your adaptive Runge-Kutta solver using again the circular orbit of the Earth from Part a). Part a) you expect that the adaptive scheme gives you the ‘same’ results. Document your results. Using the same initial conditions and the same criteria that you had used in After feeling c Use as an initial test the following parameters initial tangential velocity isomfortable with your program you will now use it for elliptical orbits. π/2 AU/yr. The initial time step is: The initial radial distance is 1 AU, and the Δt = 0.1 yr. Use the f ‘safety factors’ S1 = 0.9 and S2 = 4.0, max attempts to reach ideal error = 100ollowing parameters: ideal fractional error Run your program for at least one full orbit and plot the adaptive timestep ΔI = 1. E-3 (this is 0.001 in scientific notation) (in years) as a function of the distance from the Sun (in AU). Next, we want to use our adaptive Runge-Kutta program to determine the orbits of actual comets. For these orbits, we need to change our calculation of the initial conditions a little bit. Let’s look at the total energy of the planet (comet)
! For a circular orbit, we have seen that^ E^ =^12 mv^2 "^ GM^ rsm
or^ mv^ r^^2 =^ GM^ r^2 S^ m ! With this, the total energy in a circular orbit is^ v^ =^ GMs^ / r ! In an elliptical orbit, the semimajor and semiminor axis, a and b, are unequal (Figure 2). The^ E^ =^ "^ GM^2 rsm eccentricity is defined as !^ e^ =^1 "^ b and the distance from the sun at perihelion (closest approach) is^2 / a^2 ! The total energy of an elliptical orbit is similar to the one for a circular orbit if the radius r is^ q^ =^ (^1 "^ e ) a replaced with the semimajor axis a ! With this the orbital velocity can be expressed as^ E^ =^ "^ GM^2 asm
! In table 1, you find the orbital data for^ v^ =^ GM some selected comets. s^ #^ $^ %^2^ r^ "^1^ a &^ '^ (^ =^4 )^2 #^ $^ %^2^ r^ "^1^ a &^ '^ ( For the comets in the table use your program to determine the time that it takes the comet to orbit once around the sun (the length of a comet year). Compare your results with those given in the table.
Figure 2. Elliptical orbit around the Sun.
r
x
y
2a
q
2b Q
Problem 4) We know from our daily weather forecast The Lorenz Model that weather is intrinsically unpredictable. This is not only because of the complexity of the system, but because of the nonlinear nature of the governing equations. Ed Lorenz, an MIT meteorologist, formulated a simple model of the global weather, which was based on a 12-variable system of nonlinear ODEs. He observed that the solution was extremely sensitive to the initial conditions. This sensitivity became known as the ‘Butterfly Effect’. We can study this effect a little bit more easily using a simpler model with only three variables. With three variables the Lorenz Model is !
dx dt = "( y # x ) !
dy dt = rx " y " xz ! where σ, r, and b are positive constants. Originally, these equations were developed as a model for^ dz^ dt^ =^ xy^ "^ bz buoyant c overturn, and y and z measure the horizontal and vertical temperature gradients. The parameters and b depend on the fluid properties and the geometry of the containeronvection in a fluid. Without going into much detail, x measures the rate of convective. The parameter r is σ proportional to the applied temperature gradient. Typical values are σ = 10, b = 8/3, and r = 28. Use your adaptive Runge changes in the initial conditions lead to large ch-Kutta program to solve the Lorenz Model and observe how smallanges in the solution. To do this use the above For the initial conditions use (x, y, z) = (1, 1, 20) and (x, y, z) = (1, 1, 20.01). Use a fractional accuracy of 1. E-mentioned values for the parameters-3 and an initial timestep of Δ σt =1. = 10, b = 8/3, and r = 28. Plot the solution of x(t) for the first 20 seconds.
