Numerical ODE Solution with Adams-Bashforth-Moulton - Prof. Donald G. Luttermoser, Study notes of Physics

Download Numerical ODE Solution with Adams-Bashforth-Moulton - Prof. Donald G. Luttermoser and more Study notes Physics in PDF only on Docsity! PHYS-4007/5007: Computational Physics Course Lecture Notes Section IX Dr. Donald G. Luttermoser East Tennessee State University Version 4.1 Abstract These class notes are designed for use of the instructor and students of the course PHYS-4007/5007: Computational Physics taught by Dr. Donald Luttermoser at East Tennessee State University. Donald G. Luttermoser, ETSU IX–3 d) Occasionally, other letters will be used instead of x for the dependent variable and y for the dependent variable (i.e., function) — the meaning will be clear from the context. e) We shall assume that it is always possible to solve a given ODE for the highest derivative, obtaining y(n) = f ( x, y, y′, y′′, ..., y(n−1) ) . (IX-6) f) A second important classification of ODEs is according to whether they are linear or nonlinear. The DE in Eq. (IX-5) is said to be linear if F is a linear function of the variables y, y′, ..., y(n). Thus the general linear ODE of order n is a0(x)y (n) + a1(x)y (n−1) + · · ·+ an(x)y = g(x). (IX-7) g) An equation which is not of the form in Eq. (IV-7) is a nonlinear equation. For example, the equation for the motion of a pendulum is nonlinear: d2θ dt2 + g ` sin θ = 0. (IX-8) 3. If a DE depends upon several independent variables, it is called a partial differential equation (PDE) =⇒ the DE contains par- tial derivatives (e.g., ∂/∂t) instead of ordinary derivatives (e.g., d/dt). The wave equation is a good example of a PDE: a2 ∂2u(x, t) ∂x2 = ∂2u(x, t) ∂t2 . (IX-9) B. Numerical Methods 1. Terminology: a) Shooting or marching methods: The solution is cal- culated step by step by starting at one boundary and in- tegrating toward the other. The step size is the change IX–4 PHYS-4007/5007: Computational Physics (e.g., ∆x) in independent variable used in a shooting or marching scheme. b) Iterative method: A repetitive process by which suc- cessively more accurate approximations to the solution are obtained. An iteration is one cycle of the repetitive process. c) Difference equation: An approximation to a DE where a derivative is replaced by a quotient. d) Relaxation method: A solution is calculated every- where at once by solving a set of difference equations in an iterative fashion. e) Computational mesh or grid: The independent vari- able is represented by a set of discrete values (e.g., a set layers of given thickness in a stellar atmosphere) called grid points, zones, or cells. The ∆xi (or ∆ri, ∆Mi, etc.) is called the grid spacing or mesh size or inter- val. ∆xi ≡ xi+1 − xi. (IX-10) f) A model then becomes the set of physical properties spec- ified at all the grid points, e.g., {Pi, Ti, Fi, ρi} at zones xi. (IX-11) g) An evolution is a sequence of models at different times tn. Each model is a generation or cycle. Each successive advance forward in time is called the time step and ∆t ≡ tn+1 − tn (IX-12) is called the time step size. Donald G. Luttermoser, ETSU IX–5 h) Truncation error (TE) is the per step (for shooting schemes) or per mesh (for relaxation schemes) error in the calculation of the dependent variables. Cumulative truncation error (CE) is the total such error across the grid. i) Roundoff error is error introduced by the finite number of digits carried by the computer. The higher precision you use, the smaller the roundoff error. j) The order of a numerical scheme is the power of the mesh size or step size in the highest order terms which are ac- curately represented by a numerical scheme. Some books use the TE, some the CE, to define this, and so there is often some confusion. k) Explicit schemes: One where values at the next step are obtained by a direct algebraic computation involving only values from the previous step. l) Implicit schemes: One where new values must be solved for interatively. 2. Expansions: a) Often there are times (especially at boundaries) where singularities appear in an equation (i.e., divide by zero). During these times, one must use expansions to avoid these formal singularities. For instance, in the hydrostatic equilibrium equation: dP dr = −ρGMr r2 , (IX-13) as r → 0, Mr → 0, so the RHS → 0/0! IX–8 PHYS-4007/5007: Computational Physics where f (3) is the 3rd derivative of f(x) and x−h ≤ ζ ≤ x + h. iii) The TE is now quadratic in h, a big improvement over Eq. (IX-17). iv) From this formalism, it can be shown that the Taylor expansion of the second derivative is f ′′(x) = f(x + h) + f(x− h)− 2f(x) h2 − 1 12 h2f (4)(ζ), (IX-21) where x−h ≤ ζ ≤ x+h. Again, the TE is quadratic in h. e) Selecting values of h: What do you pick for h? i) First, define the absolute error as ε = |(true value)− (computed value)|. (IX-22) ii) Neglecting roundoff error, to make the TE term in Eq. (IX-20) small with respect to the other term in this equation, then choose h < √√√√ 6ε |f (3)(ζ)| . (IX-23) iii) Generally, we don’t know f (3), but often we can set a bound. For example, if f(x) = sin(x), then |f (3)(ζ)| ≤ 1 so if we want an absolute error of ε ≈ 10−6 (a typical value), then we should take h ≈ 2× 10−3. iv) If we cannot estimate an upper bound, then arbi- trarily pick a value of h, use it, try a smaller value Donald G. Luttermoser, ETSU IX–9 of h, compare the two answers, and if they are close enough, assume everything is fine and your answer is converged. If not, keep on choosing smaller val- ues of h (i.e., iterate) until the above is true. This is not universally true however! 3. Shooting Methods: a) Assume we have a DE given by dy/dx = f(x, y) as shown in Figure (IX-1). Now, given (xj, yj) for all j ≤ i, obtain (xi+1, yi+1). This is the standard technique in the shoot- ing or marching method. The are a variety of ways to carry out such a calculation. b) Euler’s Method: Assume we wish to follow the motion of a mass m using Newton’s 2nd Law of Motion. The equation of motion that we want to solve numerically is d~v dt = ~a(~r,~v) = 1 m ~F, (IX-24) where d~r dt = ~v (IX-25) and ~a is the acceleration. Euler’s method uses the right derivative formula (see Eq. IX-18), where we replace the grid step h with the time step τ . The equations of motion become ~v(t + τ)− ~v(t) τ + O(τ) = ~a[~r(t), ~v(t)] (IX-26) ~r(t + τ)− ~r(t) τ + O(τ) = ~v(t) (IX-27) or ~v(t + τ) = ~v(t) + τ~a[~r(t), ~v(t)] + O(τ 2) (IX-28) ~r(t + τ) = ~r(t) + τ~v(t) + O(τ 2). (IX-29) IX–10 PHYS-4007/5007: Computational Physics Figure IX–1: Data points given by the DE dy/dx = f(x, y). Notice that τO(τ) = O(τ 2). i) We introduce the notation: fn ≡ f [(n− 1)τ ]; n = 1, 2, ... (IX-30) so f1 = f(t = 0). Our equations for the Euler method (dropping the error term) now take the form ~vn+1 = ~vn + τ~an (IX-31) ~rn+1 = ~rn + τ~vn. (IX-32) ii) The calculation of trajectory would proceed as follows: 1. Specify the initial conditions ~ri and ~vi. 2. Choose a time step τ . 3. Calculate the acceleration given the current ~r and ~v. 4. Use the Euler method to compute the new ~r and ~v. 5. Go to step 3 until enough trajectory points have been computed. Donald G. Luttermoser, ETSU IX–13 a) Grid: First, one needs to set up a grid of independent variable values from one boundary to the other. Let’s assume that ri corresponds to this grid, then the grid is represented as ri, i = 0, 1, ..., N or N + 1 grid points. b) Model: The model is defined as the set of dependent variables at each grid point, e.g., {Pi, Ti, Mri, Lri} at ri for 4N + 4 unknowns. c) Then the Differential Equations take on the form dP dr = f1(P, T, Mr, Lr, r) (IX-42) dT dr = f2(P, T, Mr, Lr, r) (IX-43) dMr dr = f3(P, T, Mr, Lr, r) (IX-44) dLr dr = f4(P, T, Mr, Lr, r). (IX-45) d) Differencing Schemes: i) Forward differencing: Pi+1 − Pi ∆ri = f1 (Pi, Ti, · · ·) . (IX-46) ii) Centered differencing: Pi+1 − Pi ∆ri = f1 ( Pi+1 + Pi 2 , Ti+1 + Ti 2 , · · · ) . (IX-47) iii) Backward differencing: Pi+1 − Pi ∆ri = f1 (Pi+1, Ti+1, · · ·) . (IX-48) IX–14 PHYS-4007/5007: Computational Physics e) Boundary Conditions: i) Center (i.e., boundary 1): C1 (P0, T0, Mr0, Lr0) = 0 C2 (P0, T0, Mr0, Lr0) = 0. (IX-49) ii) Surface (i.e., boundary 2): S1 (PN , TN , MrN , LrN) = 0 S2 (PN , TN , MrN , LrN) = 0. (IX-50) f) Difference Equations: i) Choose one of the differencing schemes and apply to Eqs. (IX-42) through (IX-45). ii) Plug in known ri values which results in 4N dif- ference equations, which can be nonlinear algebraic equations for 4N + 4 unknowns. iii) These equations can be written formally as g1 (Pi, Ti,Mri, Lri, Pi+1, Ti+1,Mri+1, Lri+1) = 0 (IX-51) g2 (Pi, Ti,Mri, Lri, Pi+1, Ti+1,Mri+1, Lri+1) = 0 (IX-52) g3 (Pi, Ti,Mri, Lri, Pi+1, Ti+1,Mri+1, Lri+1) = 0 (IX-53) g4 (Pi, Ti,Mri, Lri, Pi+1, Ti+1,Mri+1, Lri+1) = 0 . (IX-54) iv) We have such equations for i = 0, 1, ..., N − 1. (IX-51) to (IX-54) 4N + 4 equations plus =⇒ in (IX-49) and (IX-50) 4N + 4 unknowns v) Great! The only problem is that these equations are not linear. They can’t be solved directly. A gen- eralized Newton-Raphson iterative scheme is used to overcome this difficulty. Donald G. Luttermoser, ETSU IX–15 g) Method: Guess to solution ←− ↓ Solve linearized Add correction equations for to the old corrections guess ↑ −→ Are corrections no small enough? yes ↓ Solution One iterates until the guesses converge to some prescribed accuracy determined by the size of the corrections. jth iteration { P (j) i , T (j) i , M (j) ri , L (j) ri } jth correction { δP (j) i , δT (j) i , δM (j) ri , δL (j) ri } . i) Consider one of the difference equations. In gen- eral, gm ( P (j) i , T (j) i , ..., P (j) i+1, ... ) 6= 0. ii) Say we know { P (j) i , T (j) i , ... } and want to calculate the jth correction. To generate linear equations we can solve for the jth correction, plug { P (j) i + δP (j) i , T (j) i + δT (j) i , ... } into gm, Taylor expand, keep terms up to first order only in the jth corrections, and then set it equal to zero. IX–18 PHYS-4007/5007: Computational Physics vi) The ↔ A matrix has a banded structure and is mostly empty.     2 x 4     4 x 8   m   4 x 8   . . . m   4 x 8     2 x 4     vii) To solve for the jth correction, this matrix must be inverted. Numerically, the number of operations involved is Full (4N + 4)× (4N + 4) ∼ (4N + 4) 2 Actual Banded Matrix ∼ 82(4N + 4) h) Gaussian Elimination: Many methods now exist for inverting banded matrices. In the last section of the notes, we were introduced to the Gaussian elimination scheme for solving sets of linear equations. I here describe this technique for ODEs. Donald G. Luttermoser, ETSU IX–19 i) Start at one corner and solve one block of equa- tions for some variables in terms of the others, e.g., for upper left 2 x 4 block, δP (j) 0 = aP0δM (j) r0 + bP0δL (j) r0 + cP0 (IX-57) δT (j) 0 = aT0δM (j) r0 + bT0δL (j) r0 + cT0. ii) Store the a, b, c’s. iii) Substitute Eq. (IX-57) in the next block and re- peat (i) and (ii) above, e.g., δM (j) r0 = aMr0δM (j) r1 + bMr0δL (j) r1 + cMr0 δL (j) r0 = aLr0δM (j) r1 + bLr0δL (j) r1 + cLr0. (IX-58) δP (j) 1 = aP1δM (j) r1 + bP1δL (j) r1 + cP1 δT (j) 1 = aT1δM (j) r1 + bT1δL (j) r0 + cT1 become the new set of Eqs. (IX-57) to substitute into the next block. iv) Solving the last set of equations at the other end of the matrix gives δM (j) rN (= δM (j) ? ) and δL (j) rN (= δL(j)? ). v) Now go backwards and back substitute into the Eqs. (IX-58) using the stored a, b, c’s to find all the jth corrections. IX–20 PHYS-4007/5007: Computational Physics i) Convergence Criterion: The iterative process is said to converge when the corrections are as small as desired, e.g., all |δPi| Pi , etc. < C = const. C ∼ 10−2 to 10−6. Most codes (at least 1-D) converge after∼3 to 5 iterations. C. Fourth-Order Runge-Kutta Method 1. When numerically solving ODE, one typically rewrites 2nd and higher-order DEs into a set of 1st-order DEs: a) A 2nd-order equation such as a = d2x dt2 . (IX-59) b) This can be rewritten as 2 equations: v = dx dt (IX-60) a = dv dt . (IX-61) 2. As shown above, Euler’s method is a first-order method which is graphically represented in Figure (IX-2). 3. The Runge-Kutta method is essentially a modified Euler’s method. a) Use the derivative at one step to extrapolate the midpoint value — use the midpoint derivative to extrapolate the function at the next step (see Figure IX-3). b) Evaluates the derivative function twice at each step τ . Cu- mulative error is of order O(τ 2), a second-order method. Donald G. Luttermoser, ETSU IX–23 Figure IX–4: Fourth-order RK with the ~Fi’s in Eqs. (IX-63) to (IX-66) being represented with ∆i’s in this figure. Also ~f is given by X in this figure. 8. You may wonder, Why 4th-order and not 8th- or 23rd-order Runge-Kutta? Well, the higher order methods have better trun- cation error but also require more computation, hence, more roundoff error. The optimum, for RK schemes, is 4th order. 9. Sometimes accuracy can be improved through use of an adaptive step size. Adaptive methods are fairly easy to incorporate and I refer you to Numerical Recipes for a description on how to incorporate them. 10. The following program is taken from Numerical Recipes and is a Fortran 77 version of the 4th-order RK scheme: IX–24 PHYS-4007/5007: Computational Physics * * <<<<<<<<<<<<<<<<<<<<<<<<<<<<< RK4 - RK4 - RK4 >>>>>>>>>>>>>>>>>>>>>>>>>>> * SUBROUTINE RK4(Y, DYDX, N, X, H, YOUT, DERIVS) * * Given values for N variables Y and their derivatives DYDX known at X, use the * 4th-order Runge-Kutta method to advance the solution over an interval H and * return the incremented variables as YOUT, which need not be a distinct * array from Y. The user supplies the subroutine DERIVS(X, Y, DYDX) which * returns the derivatives DYDX at X. * IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL DERIVS PARAMETER (NMAX = 10) DIMENSION Y(N), DYDX(N), YOUT(N), YT(NMAX), DYT(NMAX), DYM(NMAX) * HH = H * 0.5D0 H6 = H / 6.D0 XH = X + HH * * First step. * DO 11 I = 1, N YT(I) = Y(I) + HH*DYDX(I) 11 CONTINUE * * Second step. * CALL DERIVS(XH, YT, DYT) DO 12 I = 1, N YT(I) = Y(I) + HH*DYT(I) 12 CONTINUE * * Third step. * CALL DERIVS(XH, YT, DYM) DO 13 I = 1, N YT(I) = Y(I) + H*DYM(I) DYM(I) = DYT(I) + DYM(I) 13 CONTINUE * * Fourth step. * CALL DERIVS(X+H, YT, DYT) DO 14 I = 1, N YOUT(I) = Y(I) + H6*(DYDX(I)+DYT(I)+2.D0*DYM(I)) 14 CONTINUE Donald G. Luttermoser, ETSU IX–25 * RETURN END D. The Adams Method: The Shampine-Gordon Routine. 1. In 1975, L.F. Shampine and M.K. Gordon wrote a book enti- tled Computer Solution of Ordinary Differential Equations: The Initial Value Problem. a) This textbook described numerical techniques in solving non-stiff initial value problems in ordinary differential equations. b) The ODE code itself is comprised of a few subroutines and a driver program as shown in the example program below. The subroutines are i) DE: Integrates a system of up to 20 first-order ODEs of the form DY(I)/DT = F(T, Y(1), Y(2), ..., Y(NEQN)) Y(I) given at T. This subroutine integrates from T to TOUT. On re- turn, the parameters in the call list are initialized for continuing the integration. The user has only to define a new value of TOUT and call DE again. ii) STEP: Integrates a system of first order ODEs over one step, normally from T to T + H, using a modified divided difference form of the Adams Pece formulas. Local extrapolation is used to im- prove absolute stability and accuracy. The code ad- justs its order at step size to control the local error per unit step in a generalized sense. Special devices are included to control roundoff error and to detect when the user is requesting too much accuracy (see program ode.f below for further details). IX–28 PHYS-4007/5007: Computational Physics c code on output so must be variables in the calling program. c c output from ode -- c c neqn -- unchanged c y(*) -- solution at t c t -- last point reached in integration. normal return has c t = tout . c tout -- unchanged c relerr,abserr -- normal return has tolerances unchanged. iflag=3 c signals tolerances increased c iflag = 2 -- normal return. integration reached tout c = 3 -- integration did not reach tout because error c tolerances too small. relerr , abserr increased c appropriately for continuing c = 4 -- integration did not reach tout because more than c 500 steps needed c = 5 -- integration did not reach tout because equations c appear to be stiff c = 6 -- invalid input parameters (fatal error) c the value of iflag is returned negative when the input c value is negative and the integration does not reach tout , c i.e., -3, -4, -5. c work(*),iwork(*) -- information generally of no interest to the c user but necessary for subsequent calls. c c subsequent calls to ode -- c c subroutine ode returns with all information needed to continue c the integration. if the integration reached tout , the user need c only define a new tout and call again. if the integration did not c reach tout and the user wants to continue, he just calls again. c the output value of iflag is the appropriate input value for c subsequent calls. the only situation in which it should be altered c is to stop the integration internally at the new tout , i.e., c change output iflag=2 to input iflag=-2 . error tolerances may c be changed by the user before continuing. all other parameters must c remain unchanged. c c*********************************************************************** c* subroutines de and step contain machine dependent constants. * c* be sure they are set before using ode . * c*********************************************************************** c logical start,phase1,nornd dimension y(neqn),work(1),iwork(5) external f Donald G. Luttermoser, ETSU IX–29 data ialpha,ibeta,isig,iv,iw,ig,iphase,ipsi,ix,ih,ihold,istart, 1 itold,idelsn/1,13,25,38,50,62,75,76,88,89,90,91,92,93/ iyy = 100 iwt = iyy + neqn ip = iwt + neqn iyp = ip + neqn iypout = iyp + neqn iphi = iypout + neqn if(iabs(iflag) .eq. 1) go to 1 start = work(istart) .gt. 0.0d0 phase1 = work(iphase) .gt. 0.0d0 nornd = iwork(2) .ne. -1 1 call de(f,neqn,y,t,tout,relerr,abserr,iflag,work(iyy), 1 work(iwt),work(ip),work(iyp),work(iypout),work(iphi), 2 work(ialpha),work(ibeta),work(isig),work(iv),work(iw),work(ig), 3 phase1,work(ipsi),work(ix),work(ih),work(ihold),start, 4 work(itold),work(idelsn),iwork(1),nornd,iwork(3),iwork(4), 5 iwork(5)) work(istart) = -1.0d0 if(start) work(istart) = 1.0d0 work(iphase) = -1.0d0 if(phase1) work(iphase) = 1.0d0 iwork(2) = -1 if(nornd) iwork(2) = 1 return end c c c subroutine de(f,neqn,y,t,tout,relerr,abserr,iflag, 1 yy,wt,p,yp,ypout,phi,alpha,beta,sig,v,w,g,phase1,psi,x,h,hold, 2 start,told,delsgn,ns,nornd,k,kold,isnold) implicit real*8(a-h,o-z) c c ode merely allocates storage for de to relieve the user of the c inconvenience of a long call list. consequently de is used as c described in the comments for ode . c c this code is completely explained and documented in the text, c computer solution of ordinary differential equations: the initial c value problem by l. f. shampine and m. k. gordon. c logical stiff,crash,start,phase1,nornd dimension y(neqn),yy(neqn),wt(neqn),phi(neqn,16),p(neqn),yp(neqn), 1 ypout(neqn),psi(12),alpha(12),beta(12),sig(13),v(12),w(12),g(13) external f c IX–30 PHYS-4007/5007: Computational Physics c*********************************************************************** c* the only machine dependent constant is based on the machine unit * c* roundoff error u which is the smallest positive number such that * c* 1.0+u .gt. 1.0 . u must be calculated and fouru=4.0*u inserted * c* in the following data statement before using de . the routine * c* machin calculates u . fouru and twou=2.0*u must also be * c* inserted in subroutine step before calling de . * c data fouru/.888d-15/ *** c*********************************************************************** c c the constant maxnum is the maximum number of steps allowed in one c call to de . the user may change this limit by altering the c following statement data maxnum/500/ c c *** *** *** c test for improper parameters c fouru = 4.0 * d1mach(4) *** if(neqn .lt. 1) go to 10 if(t .eq. tout) go to 10 if(relerr .lt. 0.0d0 .or. abserr .lt. 0.0d0) go to 10 eps = dmax1(relerr,abserr) if(eps .le. 0.0d0) go to 10 if(iflag .eq. 0) go to 10 isn = isign(1,iflag) iflag = iabs(iflag) if(iflag .eq. 1) go to 20 if(t .ne. told) go to 10 if(iflag .ge. 2 .and. iflag .le. 5) go to 20 10 iflag = 6 return c c on each call set interval of integration and counter for number of c steps. adjust input error tolerances to define weight vector for c subroutine step c 20 del = tout - t absdel = dabs(del) tend = t + 10.0d0*del if(isn .lt. 0) tend = tout nostep = 0 kle4 = 0 stiff = .false. releps = relerr/eps abseps = abserr/eps if(iflag .eq. 1) go to 30 Donald G. Luttermoser, ETSU IX–33 c per unit step in a generalized sense. special devices are included c to control roundoff error and to detect when the user is requesting c too much accuracy. c c this code is completely explained and documented in the text, c computer solution of ordinary differential equations: the initial c value problem by l. f. shampine and m. k. gordon. c c c the parameters represent: c x -- independent variable (real*8) c y(*) -- solution vector at x (real*8) c yp(*) -- derivative of solution vector at x after successful c step (real*8) c neqn -- number of equations to be integrated (integer*4) c h -- appropriate step size for next step. normally determined by c code (real*8) c eps -- local error tolerance. must be variable (real*8) c wt(*) -- vector of weights for error criterion (real*8) c start -- logical variable set .true. for first step, .false. c otherwise (logical*4) c hold -- step size used for last successful step (real*8) c k -- appropriate order for next step (determined by code) c kold -- order used for last successful step c crash -- logical variable set .true. when no step can be taken, c .false. otherwise. c the arrays phi, psi are required for the interpolation subroutine c intrp. the array p is internal to the code. all are real*8 c c input to step c c first call -- c c the user must provide storage in his driver program for all arrays c in the call list, namely c c dimension y(neqn),wt(neqn),phi(neqn,16),p(neqn),yp(neqn),psi(12) c c the user must also declare start and crash logical variables c and f an external subroutine, supply the subroutine f(x,y,yp) c to evaluate c dy(i)/dx = yp(i) = f(x,y(1),y(2),...,y(neqn)) c and initialize only the following parameters: c x -- initial value of the independent variable c y(*) -- vector of initial values of dependent variables c neqn -- number of equations to be integrated c h -- nominal step size indicating direction of integration IX–34 PHYS-4007/5007: Computational Physics c and maximum size of step. must be variable c eps -- local error tolerance per step. must be variable c wt(*) -- vector of non-zero weights for error criterion c start -- .true. c c step requires the l2 norm of the vector with components c local error(l)/wt(l) be less than eps for a successful step. the c array wt allows the user to specify an error test appropriate c for his problem. for example, c wt(l) = 1.0 specifies absolute error, c = dabs(y(l)) error relative to the most recent value of c the l-th component of the solution, c = dabs(yp(l)) error relative to the most recent value of c the l-th component of the derivative, c = dmax1(wt(l),dabs(y(l))) error relative to the largest c magnitude of l-th component obtained so far, c = dabs(y(l))*relerr/eps + abserr/eps specifies a mixed c relative-absolute test where relerr is relative c error, abserr is absolute error and eps = c dmax1(relerr,abserr) . c c subsequent calls -- c c subroutine step is designed so that all information needed to c continue the integration, including the step size h and the order c k , is returned with each step. with the exception of the step c size, the error tolerance, and the weights, none of the parameters c should be altered. the array wt must be updated after each step c to maintain relative error tests like those above. normally the c integration is continued just beyond the desired endpoint and the c solution interpolated there with subroutine intrp . if it is c impossible to integrate beyond the endpoint, the step size may be c reduced to hit the endpoint since the code will not take a step c larger than the h input. changing the direction of integration, c i.e., the sign of h , requires the user set start = .true. before c calling step again. this is the only situation in which start c should be altered. c c output from step c c successful step -- c c the subroutine returns after each successful step with start and c crash set .false. . x represents the independent variable c advanced one step of length hold from its value on input and y c the solution vector at the new value of x . all other parameters c represent information corresponding to the new x needed to Donald G. Luttermoser, ETSU IX–35 c continue the integration. c c unsuccessful step -- c c when the error tolerance is too small for the machine precision, c the subroutine returns without taking a step and crash = .true. . c an appropriate step size and error tolerance for continuing are c estimated and all other information is restored as upon input c before returning. to continue with the larger tolerance, the user c just calls the code again. a restart is neither required nor c desirable. c logical start,crash,phase1,nornd dimension y(neqn),wt(neqn),phi(neqn,16),p(neqn),yp(neqn),psi(12) dimension alpha(12),beta(12),sig(13),w(12),v(12),g(13), 1 gstr(13),two(13) external f c*********************************************************************** c* the only machine dependent constants are based on the machine unit * c* roundoff error u which is the smallest positive number such that * c* 1.0+u .gt. 1.0 . the user must calculate u and insert * c* twou=2.0*u and fouru=4.0*u in the data statement before calling * c* the code. the routine machin calculates u . * c data twou,fouru/.444d-15,.888d-15/ *** c*********************************************************************** data two/2.0d0,4.0d0,8.0d0,16.0d0,32.0d0,64.0d0,128.0d0,256.0d0, 1 512.0d0,1024.0d0,2048.0d0,4096.0d0,8192.0d0/ data gstr/0.500d0,0.0833d0,0.0417d0,0.0264d0,0.0188d0,0.0143d0, 1 0.0114d0,0.00936d0,0.00789d0,0.00679d0,0.00592d0,0.00524d0, 2 0.00468d0/ c data g(1),g(2)/1.0,0.5/,sig(1)/1.0/ c c twou = 2.0 * d1mach(4) *** fouru = 2.0 * twou *** c *** begin block 0 *** c check if step size or error tolerance is too small for machine c precision. if first step, initialize phi array and estimate a c starting step size. c *** c c if step size is too small, determine an acceptable one c crash = .true. if(dabs(h) .ge. fouru*dabs(x)) go to 5 h = dsign(fouru*dabs(x),h) return IX–38 PHYS-4007/5007: Computational Physics temp4 = k*kp1 v(k) = 1.0d0/temp4 nsm2 = ns-2 if(nsm2 .lt. 1) go to 130 do 125 j = 1,nsm2 i = k-j 125 v(i) = v(i) - alpha(j+1)*v(i+1) c c update v(*) and set w(*) c 130 limit1 = kp1 - ns temp5 = alpha(ns) do 135 iq = 1,limit1 v(iq) = v(iq) - temp5*v(iq+1) 135 w(iq) = v(iq) g(nsp1) = w(1) c c compute the g(*) in the work vector w(*) c 140 nsp2 = ns + 2 if(kp1 .lt. nsp2) go to 199 do 150 i = nsp2,kp1 limit2 = kp2 - i temp6 = alpha(i-1) do 145 iq = 1,limit2 145 w(iq) = w(iq) - temp6*w(iq+1) 150 g(i) = w(1) 199 continue c *** end block 1 *** c c *** begin block 2 *** c predict a solution p(*), evaluate derivatives using predicted c solution, estimate local error at order k and errors at orders k, c k-1, k-2 as if constant step size were used. c *** c c change phi to phi star c if(k .lt. nsp1) go to 215 do 210 i = nsp1,k temp1 = beta(i) do 205 l = 1,neqn 205 phi(l,i) = temp1*phi(l,i) 210 continue c c predict solution and differences c Donald G. Luttermoser, ETSU IX–39 215 do 220 l = 1,neqn phi(l,kp2) = phi(l,kp1) phi(l,kp1) = 0.0d0 220 p(l) = 0.0d0 do 230 j = 1,k i = kp1 - j ip1 = i+1 temp2 = g(i) do 225 l = 1,neqn p(l) = p(l) + temp2*phi(l,i) 225 phi(l,i) = phi(l,i) + phi(l,ip1) 230 continue if(nornd) go to 240 do 235 l = 1,neqn tau = h*p(l) - phi(l,15) p(l) = y(l) + tau 235 phi(l,16) = (p(l) - y(l)) - tau go to 250 240 do 245 l = 1,neqn 245 p(l) = y(l) + h*p(l) 250 xold = x x = x + h absh = dabs(h) call f(x,p,yp) c c estimate errors at orders k,k-1,k-2 c erkm2 = 0.0d0 erkm1 = 0.0d0 erk = 0.0d0 do 265 l = 1,neqn temp3 = 1.0d0/wt(l) temp4 = yp(l) - phi(l,1) if(km2)265,260,255 255 erkm2 = erkm2 + ((phi(l,km1)+temp4)*temp3)**2 260 erkm1 = erkm1 + ((phi(l,k)+temp4)*temp3)**2 265 erk = erk + (temp4*temp3)**2 if(km2)280,275,270 270 erkm2 = absh*sig(km1)*gstr(km2)*dsqrt(erkm2) 275 erkm1 = absh*sig(k)*gstr(km1)*dsqrt(erkm1) 280 temp5 = absh*dsqrt(erk) err = temp5*(g(k)-g(kp1)) erk = temp5*sig(kp1)*gstr(k) knew = k c c test if order should be lowered c IX–40 PHYS-4007/5007: Computational Physics if(km2)299,290,285 285 if(dmax1(erkm1,erkm2) .le. erk) knew = km1 go to 299 290 if(erkm1 .le. 0.5d0*erk) knew = km1 c c test if step successful c 299 if(err .le. eps) go to 400 c *** end block 2 *** c c *** begin block 3 *** c the step is unsuccessful. restore x, phi(*,*), psi(*) . c if third consecutive failure, set order to one. if step fails more c than three times, consider an optimal step size. double error c tolerance and return if estimated step size is too small for machine c precision. c *** c c restore x, phi(*,*) and psi(*) c phase1 = .false. x = xold do 310 i = 1,k temp1 = 1.0d0/beta(i) ip1 = i+1 do 305 l = 1,neqn 305 phi(l,i) = temp1*(phi(l,i) - phi(l,ip1)) 310 continue if(k .lt. 2) go to 320 do 315 i = 2,k 315 psi(i-1) = psi(i) - h c c on third failure, set order to one. thereafter, use optimal step c size c 320 ifail = ifail + 1 temp2 = 0.5d0 if(ifail - 3) 335,330,325 325 if(p5eps .lt. 0.25d0*erk) temp2 = dsqrt(p5eps/erk) 330 knew = 1 335 h = temp2*h k = knew if(dabs(h) .ge. fouru*dabs(x)) go to 340 crash = .true. h = dsign(fouru*dabs(x),h) eps = eps + eps return Donald G. Luttermoser, ETSU IX–43 c xout by evaluating the polynomial there. information defining this c polynomial is passed from step so intrp cannot be used alone. c c this code is completely explained and documented in the text, c computer solution of ordinary differential equations: the initial c value problem by l. f. shampine and m. k. gordon. c c input to intrp -- c c all floating point variables are double precision c the user provides storage in the calling program for the arrays in c the call list dimension y(neqn),yout(neqn),ypout(neqn),phi(neqn,16),psi(12) c and defines c xout -- point at which solution is desired. c the remaining parameters are defined in step and passed to intrp c from that subroutine c c output from intrp -- c c yout(*) -- solution at xout c ypout(*) -- derivative of solution at xout c the remaining parameters are returned unaltered from their input c values. integration with step may be continued. c dimension g(13),w(13),rho(13) data g(1)/1.0d0/,rho(1)/1.0d0/ c hi = xout - x ki = kold + 1 kip1 = ki + 1 c c initialize w(*) for computing g(*) c do 5 i = 1,ki temp1 = i 5 w(i) = 1.0d0/temp1 term = 0.0d0 c c compute g(*) c do 15 j = 2,ki jm1 = j - 1 psijm1 = psi(jm1) gamma = (hi + term)/psijm1 eta = hi/psijm1 limit1 = kip1 - j IX–44 PHYS-4007/5007: Computational Physics do 10 i = 1,limit1 10 w(i) = gamma*w(i) - eta*w(i+1) g(j) = w(1) rho(j) = gamma*rho(jm1) 15 term = psijm1 c c interpolate c do 20 l = 1,neqn ypout(l) = 0.0d0 20 yout(l) = 0.0d0 do 30 j = 1,ki i = kip1 - j temp2 = g(i) temp3 = rho(i) do 25 l = 1,neqn yout(l) = yout(l) + temp2*phi(l,i) 25 ypout(l) = ypout(l) + temp3*phi(l,i) 30 continue do 35 l = 1,neqn 35 yout(l) = y(l) + hi*yout(l) return end <<<<<< Subroutine F >>>>>> SUBROUTINE F(X, Y, YP, ERRMSG) * * This subroutine sets up the ODEs for the quasi-static solar coronal loop * models. * IMPLICIT DOUBLE PRECISION (A-H, O-Z) REAL*8 KAPPA0 CHARACTER ERRMSG*78, INMSG*78 * DIMENSION Y(3), YP(3) * PARAMETER (NPTS=100) COMMON /CURRENT/ NZCURR COMMON /VARS/ FCOND(NPTS), T(NPTS), ELECN(NPTS), S(NPTS), Z(NPTS), 1 AREA(NPTS), GS(NPTS), E_INPUT(NPTS), RADLOSS(NPTS), NZONES * * Define constants. * KAPPA0 = 1.D-6 BK = 1.380662D-16 AMU = 1.6605654D-24 Donald G. Luttermoser, ETSU IX–45 HMASS = 1.008D0*AMU * * ERRMSG on input contains information concerning the subroutine that called F. * INMSG = ERRMSG ERRMSG = ’ ’ * * Check for fatal errors: T (Y(1)) and ELECN (Y(3)) should not be less than 0. * IF (Y(1) .LE. 0.D0) ERRMSG = 1 ’Temperature is less than or equal to 0!’ IF (Y(3) .LE. 0.D0) ERRMSG = 1 ’Electron density is less than or equal to 0!’ * DO 1 I = 1, 3 1 YP(I) = 0.0D0 * IF (ERRMSG(1:1) .NE. ’ ’) THEN * WRITE(6, 5) X, Y(1), Y(2), Y(3), INMSG * 5 FORMAT(’Error in F: s = ’, 1PE11.4,’, T(e) = ’,E11.4, * 1 ’, F(cond) = ’,E11.4,’,’,/12X,’N(e) = ’,E11.4,/’Called ’, * 2 ’from ’,A78) RETURN ENDIF * * Get the radiative losses. * RLOSS = RADLOOKUP(Y(1)) * * Get energy input. * EINP = ENERGY_INP(Y(1), 0) * * The temperature gradient. * YP(1) = -Y(2) / KAPPA0 / Y(1)**2.5 / AREA(NZCURR) * * The conductive flux gradient. * YP(2) = AREA(NZCURR) * (EINP - RLOSS*Y(3)**2) * * The electron density gradient. * YP(3) = (Y(3)/Y(1)) * ((HMASS*GS(NZCURR)/2.D0/BK) - YP(1)) * RETURN END