Thompson - Computing for Scientists and Engineers (523188), страница 57
Текст из файла (страница 57)
There an as small as 10-6 was found to make adrastic change in the quantum-oscillator wave function. We now revisit this problem, but for the simpler example of the differential equation (8.91) having F(x) = 1and initial values either increasing (stable behavior) or decreasing (unstable behavior).
We use the Noumerov method, for which the increasing exponential was investigated numerically in Exercise 8.26.Exercise 8.33(a) Compile Numerical DE_2; Noumerov Method so that FUNC retums(-x) . Make an input optionvalue = -1 and ANALYT returns value = expThenso that the second starting value can be multiplied by the factorrun the program for the Noumerov method with xmin = 0, xmax = 1, usingstepsize h = 0.05. Use as input values y [1] = y (xmin) = 1 and y [2] =Examine the computed y [k] and prepare a graph to compare the analytic and numerical-solution values for values in the range 10-4 to 10-2.(b) Repeat the above calculations for stepsizes of h = 0.025 and 0.01, goingout to x = 1 each time. Show that the strong sensitivity to the accuracy of thesecond starting value persists for these smaller values of the stepsize.
nThere are two remarks to be made about such an investigation of this very stiffsecond-order differential-equation solution. First, it doesn’t matter at what x valuethe effect of a small error is introduced, since exp (x) = exp (x1) exp (x - x1)for any value of x1, so only the distance the error has been propagated is significant.Second, this differential equation is particularly stiff because all the informationabout the initial slope is contained in the relative values of y [1] and y [2] .
Itwould really be better, if possible, to have other ways of indicating whether it is theincreasing or the decreasing solution that is to be found. For example, it might bepossible to provide accurate first derivatives near the initial values.There must be better ways to handle such stiff differential equations. The logarithmic transformation of the dependent variable that we now consider is such a way.The Riccati transformationConsider the logarithmic transformation of y in the differential equation (8.91), thatis, let(8.95)8.6 INTRODUCTION TO STIFF DIFFERENTIAL EQUATIONS307FIGURE 8.20 The function Y(x) in (8.96) is shown by the solid curve, the exponential function is indicated by the dashed curve, while the Riccati transformation of y, R (x), is shown by thedotted curve.The transformation (8.95) is called the Riccati transformation of y (x).
Notethat the second relation is more general than the first, because it allows negative values of y. It is therefore to be preferred, although the first relation is the one mostcommonly encountered. The second relation has its own problems because yshould not be identically zero anywhere it is used.The motivation for this transformation is shown graphically in Figure 8.20 forthe function(8.96)with a = 1 and b = 0.5. The choice of a and b in (8.96) is such that y (0) = 1 andy' (0) = -1. Note that on the linear scale the different slope of y from the exponentially decaying function is scarcely distinguishable, which is why a differentialequation solver will have trouble distinguishing between them.
On the other hand,the “logarithmic derivative” R (x) = y'/y, shown in the lower part of Figure 8.20, clearly distinguishes the exponential decay of slope -a from the modulating effect of the trigonometric factor.Exercise 8.34(a) Differentiate (8.96) with respect to x, then divide the derivative throughoutby y in order to show that(8.97)so that in the logarithmic derivative there is no exponential decay.308SECOND-ORDER DIFFERENTIAL EQUATIONS(b) Prepare graphs of the logarithmic derivative (8.97) for b = 0.25, 0.5, and0.75 with a = 1 over the range 0 < x < 4.
Compare with Figure 8.20 and justifywhether it is a or b that is the more significant parameter in determining thebehavior of the function R (x). nAfter this pictorial introduction to the Riccati transform we need some analysis toformulate the differential equation satisfied by R (x). Suppose that y is a solution of(8.98)This can be rewritten in terms of R as(8.99)Therefore, solving the second-order differential equation (8.98) can be divided intothe solution of two first-order equations. First, (8.99) is solved for R, then (8.95)is solved for y in terms of R.Exercise 8.35Derive the first-order Riccati form of the differential equation (8.99) for R (x),starting from (8.98).
nNow that we have the mathematical analysis in good shape, we should try somenumerics and programming to check out the Riccati transformation.Programming the Riccati algorithmThe program Numerical DE; Second order by Riccati is based on the program Numerical DE_2; Euler-type Methods in Section 8.4. Here is the Riccati program.PROGRAM 8.3 Solution of second-order differential equations by the Riccati transform.#include <stdio.h>#include <math.h>#define MAX 201main()/* Numerical DE; Second order by Riccati */FILE *fout;FILE *fopen();double Rb[MAX],errRb[MAX],y[MAX],erry[MAX];double xmin,xmax,h,twoh,a,b,x,analRb,analy;int choice,Nk,nmax,k;8.6INTRODUCTION TO STIFF DIFFERENTIAL EQUATIONSchar wa;double FUNCRb(),FUNCy(),ANALYTRb(),ANALYTy();printf("Numerical DE; Second Order by Riccati\n");printf("Write over output (w) or Add on (a): ");fout = fopen("NUMDE.NumRic.l",&wa);scanf("%s",&wa) ;printf("Input xmin, xmax, h: ");scanf("%le%le%le",&xmin,&xmax,&h) ;Nk = (xmax-xmin)/h+l.l;if(Nk > MAX-1){printf("!! # of xsteps, %i, > %i\n",Nk,MAX-1); exit(l);}twoh = 2*h;printf("Input a,b: ") ;scanf("%le%le",&a,&b) ;Rb[l] = a+b; Rb[2] = (a+b-h)/(1+(a+b)*h);x=xmin;/* Euler central for Rb */for(k=2;k<Nk;k=k+2)x = x+h;Rb[k+l] = Rb[k-l]+twoh*FUNCRb(Rb,k);Rb[k+2] = Rb[k]+twoh*FUNCRb(Rb,k+l);x = x+h;}Y[l] = 1; y[2] = 1+b*h;x = xmin;/* Euler central for y */for(k=2;k<Nk;k=k+2)x = x+h;y[k+1] = y[k-1]+twoh*FUNCy(a,Rb,y,k);y[k+2] = y[k]+twoh*FUNCy(a,Rb,y,k+l);x = x+h;}printf("Comparison with analytical values\n");x = xmin;for ( k = 1; k <= Nk; k++ )analRb = ANALYTRb(x,a,b);/* Rb analytical & error */errRb[k] = analRb-Rb[k];analy = ANALYTy(x,a,b);/* y analytical & error */erry[K] = analy-y[k] ;printf("%5.2le %6.3lg %6.3lg %6.3lg %6.3lg %6.3lg\n",x, analRb, errRb[k] , analy,y [k] , erry [k] ) ;fprintf(fout,"%5.2le %6.3lg %6.3lg %6.3lg %6.3lg %6.3lg\n",309310SECOND-ORDER DIFFERENTIAL EQUATIONSx,analRb,errRb[k],analy,y[k],erry[k]);x = x+h;printf("\nEnd Numerical DE; Second Order by Riccati");double FUNCRb(Rb,k)/* dRb/dx = FUNCRb */double Rb[];int k;double value;value = Rb[k]*Rb[k]-1;return value;double FUNCy(a,Rb,y,k)/* dy/dx = FUNCy */double Rb[],y[];double a;int k;{double value;value = (Rb[k]-a)*y[k] ;return value;}double ANALYTRb(x,a,b)/* Analytical solution for Rb=y'/y+a for oscillator */double x,a,b;double value;value = (-sin(x)+(a+b)*cos(x))/(cos(x)+(a+b)*sin(x));return value;double ANALYTy(x,a,b)/* Analytical solution for damped oscillator */double x,a,b;{double value;value = exp(-a*x)*(cos(x)+(a+b)*sin(x));return value;}8.6 INTRODUCTION TO STIFF DIFFERENTIAL EQUATIONS311The structure of Numerical DE; Second order by Riccati is straightforward and the coding is very similar to that in the Euler-type and Noumerov methods,so you can proceed directly to implement, test, and use the program.Exercise 8.36(a) Codetheprogram Numerical DE; Second order by Riccati to compute numerically the solution of (8.98).
Test the program with F(x) = F (aconstant).(b) Use your program to compute the numerical solution of the damped-exponential equation that has solution (8.96). Compare your numerical results fora = 1 and b = 0.5 with the analytical solution (8.96) and with the graphicaloutput in Figure 8.20. nWe now have some idea how stiff differential equations of the exponentially decreasing type can be solved numerically, so it is interesting to consider brieflysecond-order differential equations with solutions that oscillate rapidly.Madelung’s transformation for stiff equationsSuppose that in the linear second-order differential equation(8.100)we have F (x) << 0 for a wide range of x. The solutions of such a differentialequation will be rapidly oscillating and therefore numerically troublesome, just aswhen F has the opposite sign and is large in magnitude we get a rapidly changingstiff differential equation that can be handled numerically by the Riccati transformation that we developed in the two preceding subsections.The Madelung transformation is designed to reduce such problems by dividing yinto an amplitude part, r, and a phase part, according to(8.101)where r is purely real.
In order to identify real and imaginary quantities, we write(8.102)then insert this and (8.101) in the differential equation (8.100). The real quantities randthen satisfy(8.103)and(8.104)312SECOND-ORDER DIFFERENTIAL EQUATIONSThese equations can be simplified in the common case that F is purely real, asyou may wish to prove.Exercise 8.37(a) Verify equations (8.103) and (8.104) by making the indicated substitutions.(b) Show that if F is purely real, then (8.104) can be integrated to give(8.105)where B is a constant of integration.(c) By substitution of this result into (8.103), show that, again only for real F,(8.106)so that if this equation can be integrated the resulted can be substituted into thefirst-order differential equation (8.105) for(d) Verify that if F (x) =with real, then r is a constant and, as expected,nTherefore the Madelung transformation for real and negative F in (8.100) can remove the stiffness of a linear differential-equation solution in much the same way asthe Riccati transformation does for large real-positive values of F.Our introduction to the mathematical and numerical analysis of differential equations should set you on the path to reading and understanding the extensive technicalliterature.
A compendium of formulas is provided in the books by Jain and by Zill,while many formulas are listed in Abramowitz and Stegun. The problem of stiff differential equations is given particular attention in the book by Gear.REFERENCES ON SECOND-ORDER EQUATIONSAbramowitz, M., and I. A. Stegun, Handbook of Mathematical Functions, Dover,New York, 1964.Backus, J., The Acoustical Foundations of Music, Norton, New York, secondedition, 1977.Bernoulli, J., 0pera Omnia, Marc-Michel Bousquet, Laussanne, 1742, reprinted byGeorge Olms, Hildesheim, Germany, 1968, edited by J. E.
Hofmann; Vol. IIpp. 232, 251, Vol. IV p. 234.Brandt, S., and H. D. Dahmen, Quantum Mechanics on the Personal Computer,Springer-Verlag, Berlin, 1989.Brandt, S., and H. D. Dahmen, Quantum Mechanics on the Macintosh, SpringerVerlag, Berlin, 199 1.REFERENCES ON SECOND-ORDER EQUATIONS313Braun, M., C. S. Coleman, and D. A. Drew (eds.), Differential Equation Models,Springer-Verlag, New York, 1983.Cohen-Tannoudji, C., B. Diu, and F.
Laloë, Quantum Mechanics, WileyInterscience, New York, 1977.Galilei, G., Two New Sciences, Elzevirs, Leyden, The Netherlands, 1638,translated by Stillman Drake, University of Wisconsin Press, Madison, 1974,p. 143.Gear, C. W., Applications and Algorithms in Computer- Science, Science ResearchAssociates, Chicago, 1978.Gilbert, D.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
