Nash - Scientific Computing with PCs (523165), страница 47
Текст из файла (страница 47)
Dobb’s JournalIf we assume that the 4th coefficient in the MINITAB solution is zero and regularize, as with MATLAB, sothat the first height is zero, we obtain heights at 0, 99.995, 121.015, and 96.02 meters (as we should!). Wehave simply subtracted -96.02 from each coefficient.Finally, we may choose or be forced to solve the least squares problem iteratively by a techniquespecifically designed to minimize the sum of squares function. These techniques are finding growingimportance in problems in statistics or geodesy where very large numbers of parameters b (into severalthousands) must be estimated. Figure 15.5.8 shows the results of applying Nash J C (1990d), Algorithm24 to this problem.
Adding 79.2575 to all elements of the solution gives benchmark heights of 0, 99.995,121.015, and 96.02 meters as anticipated.Figure 15.5.6MTBMTBMTBMTBMTB>>>>>MTB >MTB >MTB >MTB >MTB >SUBC>**Partial MINITAB output for the height difference problem 15.1.1bexec ’survey.xec’note MINITAB attempt to carry outnote Nash J C CNM2 Surveying example pg 240note matrix defines which height differences are involvedread c1-c45 ROWS READend datanote the measured height differencesset c5 (data is entered here)end dataregr c5 on 4 c1-c4;noconstant.C4 is highly correlated with other X variablesC4 has been removed from the equationThe regression equation isC5 = - 96.0 C1 + 3.97 C2 + 25.0 C3PredictorCoefNoconstantC1-96.0200C23.97500C324.9950s = 0.01732Figure 15.5.7Stdevt-ratiop0.01730.013690.0137-5543.65290.291825.350.0000.0000.000Edited MATLAB output for the height difference problem.
The names of matrices andvectors differ from those used in the text.% Nash CNM2 Surveying example pg 240% SURVEY.M% matrix defines which height differences areinvolvedA = [ 1 -1 0 0;0 1 -1 0;0 0 1 -1;1 0 -1 0;0 1 0 -1]% and the measured height differencesy = [ -99.99; -21.03; 24.98; -121.02; 3.99]% true valuesbtrue = [0; 100; 121; 96]% compute current residuals from truedifferencesw=A*btruer=y-w = 0.0100 -0.0300 -0.0200 -0.0200-0.0100SS=r’*r = 0.0019%initial estimates from first 3 data pointsAA=A(1:3,1:3)yy=y(1:3)xx=AA\yyxp=[xx;0]xp=xp-ones(4,1)*xp(1)xp = 0 99.9900 121.0200 96.0400%errorsdd=xp-btruedd = 0 -0.0100 0.0200 0.0400%sum of squaresdd’*dd = 0.0021%residuals from initial estimatesr=y-A*xp = 0 0 -0.0000 0 0.0400ss=r’*r = 0.0016% solve, noting singularityx = A\yWarning: Rank deficient, rank = 3 tol =1.9230e-015x = -96.0200 3.9750 24.9950 0w=A*x;% residualsr=y-w = 0.0050 -0.0100 -0.0150 -0.0050 0.0150% and sum of squaresSS=r’*r = 6.0000e-004% regularize solution[n,k]=size(x);z=x-x(1)*ones(n,1)=0 99.9950 121.0150 96.0200w=A*z=-99.995 -21.02 24.995 -121.015 3.975% residualsr=y-w = 0.005 -0.01 -0.015 -0.005 0.015% and sum of squaresSS=r’*r = 6.0000e-004%errorsdd=z-btrue = 0 -0.0050 0.0150 0.0200%sum of squaresdd’*dd = 6.5000e-00415: DATA FITTING AND MODELLINGFigure 15.5.8137Partial output of iterative method for height difference data, part (b) of Table 15.1.1.dr24ls.pas -- linear least squares by conjugate gradients1992/08/11 18:24:36Number of rows in coefficient matrix = 5Number of columns in coefficient matrix = 4Matrixin.pas -- generate or input a real matrix 5 by 4Row 11.00000-1.000000.000000.00000Row 20.000001.00000-1.000000.00000Row 30.000000.000001.00000-1.00000Row 41.000000.00000-1.000000.00000Row 50.000001.000000.00000-1.00000RHS vector-9.9990E+01 -2.1030E+01 2.4980E+01 -1.2102E+02 3.9900E+00Initial guess for solution -- all zerosNormal equations -- coefficient matrix2.00000-1.00000-1.000000.00000-1.000003.00000-1.00000-1.00000-1.00000-1.000003.00000-1.000000.00000-1.00000-1.000002.00000Normal equations - RHS-221.0100082.95000 167.03000 -28.97000Step too large -- coefficient matrix indefinite?Solution after 4 iterations.
Est. normal eqn. sumsquares -1.0000000000E+35-79.2575020.7375041.7575016.76250For original least squares problem -- Residuals-5.000E-03 1.000E-02 1.500E-02 5.000E-03 -1.500E-02Sum of squared residuals = 6.0000000025E-0415.6 AssessmentFor the growth curve problem, the choice among the four general styles of solution tool, that is, statisticalpackage, spreadsheet add-in, special purpose program or package or general mathematical problemsolving environment, will depend on the user’s familiarity with the tools, the expected number of timessimilar problems need to be solved, the additional output required, and the degree of difficulty theproblem presents.
For the weed data, bad scaling and the singularity of the Hessian means that usuallygood techniques, such as Newton’s method, may fail or will appear inefficient. Supposedly inferiormethods can make progress from points where a straightforward Newton method cannot go on. It ispossible, of course, to adjust the Newton method to overcome such difficulties. However, the user mustnow cope with even more detail in the description of the method, without assistance of off-the-shelfsoftware.Previous Home Next138Copyright © 1984, 1994 J C & M M NashNash Information Services Inc., 1975 Bel Air Drive, Ottawa, ON K2C 0X1 CanadaSCIENTIFIC COMPUTING WITH PCsCopy for:Dr.
Dobb’s JournalChapter 16Trajectories of Spacecraft16.116.216.316.416.516.6Problem statementFormulationsMethods and AlgorithmsPrograms or PackagesSome solutionsAssessmentThis chapter looks at the orbits spacecraft follow near the earth and moon. The problem is a goodillustration of differential equation solution.16.1 Problem StatementTrajectories of spacecraft are an important class of problems that can be solved by numerical methods.Consider the problem of the Apollo capsule in the U.S. space program of the late 1960s. The system ofinterest involves the earth, moon, and the spacecraft.
We will restrict our attention to trajectories in theplane of rotation of the moon about the earth to simplify our work (Figure 16.1.1).16.2 FormulationsForsythe, Malcolm, and Moler (1977, pg. 153-4) and Kahaner, Moler and Nash S G (1989, p. 333) presenta version of this problem, but leave out important details of constants and units. The simplified problemis more difficult to relate to the real world, though the equations are simpler.
Here we shall sketch aderivation of the appropriate equations based on classical mechanics (Synge and Griffith, 1959) that resultsin a set of linked first order ordinary differential equations. There are starting values that give the speedand position of the spacecraft whose trajectory we wish to calculate, so this is an initial value problem.The trajectory is not necessarily a stable orbit (periodic trajectory). If we start with a periodic orbit andperturb them, algebraic eigenproblems are derived from the equations. Such a formulation is appropriate,for example, for earth orbits of artificial satellites influenced by moon gravitation and perturbations ofplanetary orbits by other planets.The main equation is Newton’s second law of motion(16.2.1)F = w a = w d 2 rs / dt 2where F is the (vector) force acting on the spacecraft, w its mass (we shall use M and m for the massesof the earth and moon respectively).
a is the spacecraft acceleration, and rs its vector position in astationary coordinate system. Time is represented by the variable t. We are more likely to be interestedin motion of the spacecraft relative to the earth and moon, rather than some abstract but stationary (thatis, non-accelerated) coordinate frame. The desired coordinates have these two "large" bodies fixed. Inreality the earth-moon coordinate system is rotating, so we introduce corrections to the usual equation ofmotion 16.2.1. We choose as our coordinate system the rectangular axes in x, y, and z directions definedas follows. The x-axis is taken along the earth-moon line of centers.
The y-axis is taken perpendicular to16: SPACECRAFT TRAJECTORIESFigure 16.1.1139Diagram showing the variables of the spacecraft systemthe x-axis and lying in the plane of rotation of the earth and moon. The z-axis is taken to be the remainingperpendicular direction the angular velocity is a positive right-handed vector. That is, if the earth andmoon are rotating anti-clockwise in the plane of the page, z is pointing up out of the paper).The corrections to the equation of motion are called the Coriolis force and the centrifugal force. Ourcorrected equations of motion are(16.2.2)w d 2 r / dt 2= F - w (at + ac )where r is the new position vector relative to the rotating axes, and at and ac are the centrifugal andCoriolis accelerations respectively.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
