1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
|
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 3.2//EN">
<!--Converted with LaTeX2HTML 97.1 (release) (July 13th, 1997)
by Nikos Drakos (nikos@cbl.leeds.ac.uk), CBLU, University of Leeds
* revised and updated by: Marcus Hennecke, Ross Moore, Herb Swan
* with significant contributions from:
Jens Lippman, Marek Rouchal, Martin Wilck and others -->
<HTML>
<HEAD>
<TITLE>Numerical Methods</TITLE>
<META NAME="description" CONTENT="Numerical Methods">
<META NAME="keywords" CONTENT="sun67">
<META NAME="resource-type" CONTENT="document">
<META NAME="distribution" CONTENT="global">
<META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=iso_8859_1">
<LINK REL="STYLESHEET" HREF="sun67.css">
<LINK REL="previous" HREF="node226.html">
<LINK REL="up" HREF="node197.html">
<LINK REL="next" HREF="node228.html">
</HEAD>
<BODY >
<BR> <HR>
<A NAME="tex2html2721" HREF="node228.html">
<IMG WIDTH="37" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="next" SRC="next_motif.gif"></A>
<A NAME="tex2html2719" HREF="node197.html">
<IMG WIDTH="26" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="up" SRC="up_motif.gif"></A>
<A NAME="tex2html2715" HREF="node226.html">
<IMG WIDTH="63" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="previous" SRC="previous_motif.gif"></A> <A HREF="sun67.html#stardoccontents"><IMG ALIGN="BOTTOM" BORDER="0"
SRC="contents_motif.gif"></A>
<BR>
<B> Next:</B> <A NAME="tex2html2722" HREF="node228.html">SUMMARY OF CALLS</A>
<BR>
<B>Up:</B> <A NAME="tex2html2720" HREF="node197.html">EXPLANATION AND EXAMPLES</A>
<BR>
<B> Previous:</B> <A NAME="tex2html2716" HREF="node226.html">Focal-Plane Astrometry</A>
<BR> <HR> <P>
<P><!--End of Navigation Panel-->
<H2><A NAME="SECTION000521000000000000000">
Numerical Methods</A>
</H2>
SLALIB contains a small number of simple, general-purpose
numerical-methods routines. They have no specific
connection with positional astronomy but have proved useful in
applications to do with simulation and fitting.
<P>
At the heart of many simulation programs is the generation of
pseudo-random numbers, evenly distributed in a given range:
sla_RANDOM
does this. Pseudo-random normal deviates, or ``Gaussian
residuals'', are often required to simulate noise and
can be generated by means of the function
sla_GRESID.
Neither routine will pass super-sophisticated
statistical tests, but they work adequately for most
practical purposes and avoid the need to call non-standard
library routines peculiar to one sort of computer.
<P>
Applications which perform a least-squares fit using a traditional
normal-equations methods can accomplish the required matrix-inversion
by calling either
sla_SMAT
(single precision) or
sla_DMAT
(double). A generally better way to perform such fits is
to use singular value decomposition. SLALIB provides a routine
to do the decomposition itself,
sla_SVD,
and two routines to use the results:
sla_SVDSOL
generates the solution, and
sla_SVDCOV
produces the covariance matrix.
A simple demonstration of the use of the SLALIB SVD
routines is given below. It generates 500 simulated data
points and fits them to a model which has 4 unknown coefficients.
(The arrays in the example are sized to accept up to 1000
points and 20 unknowns.) The model is:
<P ALIGN="CENTER">
<I>y</I> = <I>C<SUB>1</SUB></I> +<I>C<SUB>2</SUB>x</I> +<I>C<SUB>3</SUB>sin</I><I>x</I> +<I>C<SUB>4</SUB>cos</I><I>x</I>
</P>
The test values for the four coefficients are
<IMG WIDTH="79" HEIGHT="27" ALIGN="MIDDLE" BORDER="0"
SRC="img325.gif"
ALT="$C_1\!=\!+50.0$">,<IMG WIDTH="71" HEIGHT="27" ALIGN="MIDDLE" BORDER="0"
SRC="img326.gif"
ALT="$C_2\!=\!-2.0$">,<IMG WIDTH="79" HEIGHT="27" ALIGN="MIDDLE" BORDER="0"
SRC="img327.gif"
ALT="$C_3\!=\!-10.0$"> and
<IMG WIDTH="79" HEIGHT="27" ALIGN="MIDDLE" BORDER="0"
SRC="img328.gif"
ALT="$C_4\!=\!+25.0$">.Gaussian noise, <IMG WIDTH="55" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
SRC="img329.gif"
ALT="$\sigma=5.0$">, is added to each ``observation''.
<P><PRE>
IMPLICIT NONE
* Sizes of arrays, physical and logical
INTEGER MP,NP,NC,M,N
PARAMETER (MP=1000,NP=10,NC=20,M=500,N=4)
* The unknowns we are going to solve for
DOUBLE PRECISION C1,C2,C3,C4
PARAMETER (C1=50D0,C2=-2D0,C3=-10D0,C4=25D0)
* Arrays
DOUBLE PRECISION A(MP,NP),W(NP),V(NP,NP),
: WORK(NP),B(MP),X(NP),CVM(NC,NC)
DOUBLE PRECISION VAL,BF1,BF2,BF3,BF4,SD2,D,VAR
REAL sla_GRESID
INTEGER I,J
* Fill the design matrix
DO I=1,M
* Dummy independent variable
VAL=DBLE(I)/10D0
* The basis functions
BF1=1D0
BF2=VAL
BF3=SIN(VAL)
BF4=COS(VAL)
* The observed value, including deliberate Gaussian noise
B(I)=C1*BF1+C2*BF2+C3*BF3+C4*BF4+DBLE(sla_GRESID(5.0))
* Fill one row of the design matrix
A(I,1)=BF1
A(I,2)=BF2
A(I,3)=BF3
A(I,4)=BF4
END DO
* Factorize the design matrix, solve and generate covariance matrix
CALL sla_SVD(M,N,MP,NP,A,W,V,WORK,J)
CALL sla_SVDSOL(M,N,MP,NP,B,A,W,V,WORK,X)
CALL sla_SVDCOV(N,NP,NC,W,V,WORK,CVM)
* Compute the variance
SD2=0D0
DO I=1,M
VAL=DBLE(I)/10D0
BF1=1D0
BF2=VAL
BF3=SIN(VAL)
BF4=COS(VAL)
D=B(I)-(X(1)*BF1+X(2)*BF2+X(3)*BF3+X(4)*BF4)
SD2=SD2+D*D
END DO
VAR=SD2/DBLE(M)
* Report the RMS and the solution
WRITE (*,'(1X,''RMS ='',F5.2/)') SQRT(VAR)
DO I=1,N
WRITE (*,'(1X,''C'',I1,'' ='',F7.3,'' +/-'',F6.3)')
: I,X(I),SQRT(VAR*CVM(I,I))
END DO
END
</PRE>
<P>
The program produces the following output:
<P><PRE>
RMS = 4.88
C1 = 50.192 +/- 0.439
C2 = -2.002 +/- 0.015
C3 = -9.771 +/- 0.310
C4 = 25.275 +/- 0.310
</PRE>
<P>
In this above example, essentially
identical results would be obtained if the more
commonplace normal-equations method had been used, and the large
<IMG WIDTH="72" HEIGHT="25" ALIGN="MIDDLE" BORDER="0"
SRC="img330.gif"
ALT="$1000\times20$"> array would have been avoided. However, the SVD method
comes into its own when the opportunity is taken to edit the W-matrix
(the so-called ``singular values'') in order to control
possible ill-conditioning. The procedure involves replacing with
zeroes any W-elements smaller than a nominated value, for example
0.001 times the largest W-element. Small W-elements indicate
ill-conditioning, which in the case of the normal-equations
method would produce spurious large coefficient values and
possible arithmetic overflows. Using SVD, the effect on the solution
of setting suspiciously small W-elements to zero is to restrain
the offending coefficients from moving very far. The
fact that action was taken can be reported to show the program user that
something is amiss. Furthermore, if element W(J) was set to zero,
the row numbers of the two biggest elements in the Jth column of the
V-matrix identify the pair of solution coefficients that are
dependent.
<P>
A more detailed description of SVD and its use in least-squares
problems would be out of place here, and the reader is urged
to refer to the relevant sections of the book <I>Numerical Recipes</I>
(Press <I>et al.</I>, Cambridge University Press, 1987).
<P>
The routines
sla_COMBN
and
sla_PERMUT
are useful for problems which involve combinations (different subsets)
and permutations (different orders).
Both return the next in a sequence of results, cycling through all the
possible results as the routine is called repeatedly.
<P>
<BR>
<P>
<BR> <HR>
<A NAME="tex2html2721" HREF="node228.html">
<IMG WIDTH="37" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="next" SRC="next_motif.gif"></A>
<A NAME="tex2html2719" HREF="node197.html">
<IMG WIDTH="26" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="up" SRC="up_motif.gif"></A>
<A NAME="tex2html2715" HREF="node226.html">
<IMG WIDTH="63" HEIGHT="24" ALIGN="BOTTOM" BORDER="0" ALT="previous" SRC="previous_motif.gif"></A> <A HREF="sun67.html#stardoccontents"><IMG ALIGN="BOTTOM" BORDER="0"
SRC="contents_motif.gif"></A>
<BR>
<B> Next:</B> <A NAME="tex2html2722" HREF="node228.html">SUMMARY OF CALLS</A>
<BR>
<B>Up:</B> <A NAME="tex2html2720" HREF="node197.html">EXPLANATION AND EXAMPLES</A>
<BR>
<B> Previous:</B> <A NAME="tex2html2716" HREF="node226.html">Focal-Plane Astrometry</A>
<BR> <HR> <P>
<P><!--End of Navigation Panel-->
<ADDRESS>
<I>SLALIB --- Positional Astronomy Library<BR>Starlink User Note 67<BR>P. T. Wallace<BR>12 October 1999<BR>E-mail:ptw@star.rl.ac.uk</I>
</ADDRESS>
</BODY>
</HTML>
|
