text:
ArticlePDF Available

# Tailored freeform optical surfaces

Authors:

## Abstract and Figures

Freeform optical surfaces embedded in three-dimensional space, without any symmetry, are tailored so as to solve the archetypal problem of illumination design: redistribute the radiation of a given small light source onto a given reference surface, thus achieving a desired irradiance distribution on that surface. The shape of the optical surface is found by solving a set of partial nonlinear differential equations. For most cases, a few topologically distinct solutions exist, given suitable boundary conditions.
Content may be subject to copyright.
Tailored freeform optical surfaces
Harald Ries
Optics & Energy Concepts, Paul-Gerhardt-Allee 42, 81245 Munich, Germany,
and Philipps-University, Renthof 5, 35032 Marburg, Germany
Julius Muschaweck
Optics & Energy Concepts, Paul-Gerhardt-Allee 42, 81245 Munich, Germany
Received January 24, 2001; accepted August 31, 2001; revised manuscript received September 4, 2001
Freeform optical surfaces embedded in three-dimensional space, without any symmetry, are tailored so as to
solve the archetypal problem of illumination design: redistribute the radiation of a given small light source
onto a given reference surface, thus achieving a desired irradiance distribution on that surface. The shape of
the optical surface is found by solving a set of partial nonlinear differential equations. For most cases, a few
topologically distinct solutions exist, given suitable boundary conditions. © 2002 Optical Society of America
OCIS codes: 080.2740; 350.4600; 080.2720.
1. INTRODUCTION
The general problem of optical design for illumination can
be stated as follows: Given a speciﬁed light source, de-
sign an optical system that redirects the light from the
source onto a target surface of given shape such as to pro-
duce a given desired irradiance distribution on the target.
In this formulation, for simplicity, multiple optical sur-
faces and gradient index distributions are not considered.
A large amount of work has been devoted to this classical
design problem. Often a simple system consisting of a
single reﬂective or refractive surface is desired. In this
work we also refer to single optical surfaces.
Heuristic approaches are based on the known optical
properties of conics. Imaging optical systems tradition-
ally rely almost exclusively on spherical surfaces. Even
nonimaging systems traditionally involve combinations of
conics. Recently we introduced the idea of tailoring an
optical surface, which consists of translating the desired
optical properties into a suitable differential equation,
which is then solved numerically to determine the shape
of the optical surface.15In this way smooth surfaces are
found, which in general are not conic sections. The pre-
vious work on tailoring was applicable to two-dimensional
systems, i.e., systems that possess translational or rota-
tional symmetry. In this paper we follow the same line of
thought for designing freeform three-dimensional sur-
faces. Although the extension from two to three dimen-
sions may seem a minor step, this is by no means true.
Two-dimensional tailoring leads to ordinary differential
equations. Such equations can easily be solved even by
heuristic approaches. Three-dimensional tailoring, in
contrast, leads to nonlinear partial differential equations
that are much more demanding.
Optimization approaches ﬁrst attempt a suitable pa-
rameterization of the surface with a ﬁnite and in general
low number of parameters. The value of these param-
eters are then to be determined by optimization based on
ray tracing and some objective function to measure how
well a surface characterized by certain values of the pa-
rameters performs. Optimization approaches are very
successful for many applications. However, if freeform
optical surfaces are to be designed for illumination tasks,
procedures based on optimization suffer from the fact that
adequate representation of optical surfaces tends to re-
quire a large number of parameters. Even slight varia-
tions of the surface shape can result in large variations in
the produced irradiance distribution. Ray tracing in-
volves a statistical element and therefore cannot assess
the optical performance precisely. At best it can provide
an estimate that becomes better as more rays are used.
ance distributions that are rich in detail, each evaluation
of the objective function with reasonable accuracy re-
quires a huge computational cost. For these reasons we
do not elaborate further on optimization approaches here.
In this paper we outline a constructive method of de-
sign that translates the desired optical performance into
a differential equation that is solved by numerical meth-
ods to yield the desired shape of the optical surface.
2. TAILORING IN THREE DIMENSIONS
Consider a ray from a small source located at saimed
at a point on the optical surface located at pas shown
in Fig. 1. The laws of reﬂection (refraction) govern
the direction into which this ray is redirected. In effect
these laws establish a relation between the surface-
normal unit vector N, the direction of the incoming ray
In Unit(ps) , and the direction of the outgoing ray
Out Unit(tp).
1n22nOut In1/2NOut nIn. (1)
Here Unit denotes the unit vector: Unit(v)v/
vv. Equation (1) applies for both reﬂection (n
1) and refraction (n1), where ndenotes the index
of refraction of the medium surrounding the source. To
analyze the power density carried by the radiation, we
need to consider an ensemble of adjacent rays. This
leads to the concept of a wave front; see Fig. 2.
590 J. Opt. Soc. Am. A /Vol. 19, No. 3 /March 2002 H. Ries and J. Muschaweck
0740-3232/2002/030590-0615.00 © 2002 Optical Society of America A wave front is a family of surfaces, each perpendicular to the local direction of the radiation. For example, the undistorted radiation emitted by a point source is de- scribed by concentric spheres centered at the location of the source. Wave fronts may be distorted by optical sur- faces or nonconstant refractive index and thereby have shapes different from spherical. The power density transported by the radiation along the direction of propa- gation is proportional to the absolute value of the Gauss- ian curvature of the wave front. The curvature tensor K is the derivative of the unit normal of a surface with re- spect to the tangent directions t1and t2: Ki j Ni tj i,j1, 2. (2) Since Nis a unit vector, its derivative is perpendicular to N, i.e., restricted to the tangent plane. The curvature is thus a two-dimensional tensor. The Gaussian curvature is the product of the two eigenvalues c1and c2of the cur- vature tensor deﬁned in Eq. (2). For free propagation, the curvature of the wave front changes along the ray such that the two curvature centers (also referred to as foci) remain ﬁxed. For the radiation of a point source, the two foci f1and f2coincide with the location of the source. For a general wave front, the position of the foci is calculated from the eigenvalues c1and c2of the curva- ture tensor: f1p 1 c1 N,f2p 1 c2 N. (3) Here we adopt the convention that positive curvature im- plies a convex surface, where the focus is behind. The ir- radiance of a freely propagating wave front changes in the direction of propagation according to Ref. 6: EpEp0 p0f1 兩兩 p0f2 pf1 兩兩 pf2 . (4) In general, the two foci for each point of the wave front form two (curved) surfaces in space where the optical sur- face has imaging properties. Equation (5) diverges if p coincides with one of the foci. In reality, sources are not true points and their brightness is not inﬁnite; therefore the irradiance approaches the brightness of the image of the source. This phenomenon is illustrated by the bright lines (caustics) formed on a sunny day at the bottom of a swimming pool. The two caustic surfaces are the gener- alization of the line caustic in two dimensions that is use- ful for illumination design.7The tailoring algorithm as- sumes a small source, which implies that the source brightness is much larger than the desired irradiance anywhere on the target. Therefore the caustic surfaces must never intersect the target surface. This fact further illustrates that our approach is rooted in nonimaging op- tics. To quantify the irradiance on the target surface, we ﬁrst need to derive how the curvature of a wave front is changed by an optical surface. To this end we compute the derivative of Eq. (2) with respect to displacements of the point pin the tangent plane of the optical surface. We are interested only in changes of the normal vectors conﬁned to the tangent plane of the optical surface: 1n22nOut InKo POut KOut JOut nPIn KIn JIn . (5) Here Kodenotes the curvature tensor of the optical sur- face, KIn the curvature tensor of the incoming optical beam at the intersection with the optical surface, and KOut the curvature tensor of the outgoing optical beam immediately after reﬂection (refraction) at the optical sur- face. The tensors PIn and POut are projection operators from the tangential plane of the incoming or the outgoing wave front onto the tangential plane of the optical sur- face. The tensors JIn and JOut are the Jacobians, which project the differentials from the tangent plane of the wave fronts onto the tangent plane of the optical surface: JIn i j ti In tjO,JOut i j ti Out tjO. (6) If we choose a coordinate system on each tangent plane such that the ﬁrst coordinate is sagittal, i.e., co-planar with the incoming and the outgoing ray, and the second coordinate is tangential, i.e., parallel on all surfaces, then the projection operators as well as the Jacobians are di- agonal and simple: Fig. 2. Curvature of the optical surface, illustrated here as a re- ﬂector, and the curvature of the incoming wave front combine to determine the curvature of the outgoing wave front. Fig. 1. Sketch of source location s, a point pon the optical sur- face, and points ton the target surface. Ndenotes the unit nor- mal to the optical surface. H. Ries and J. Muschaweck Vol. 19, No. 3/ March 2002 /J. Opt. Soc. Am. A 591 PIn JIn NIn 0 01 ; POut JOut NOut 0 01 . (7) In particular, for a reﬂector surface this leads to n1, Out In ⫽⫺cos(2 ) and NOut ⫽⫺NIn cos( ), where is the incidence angle, and Eq. (5) be- comes 21cos2 兲兲Ko 10 0 cos 1 KIn KOut cos 0 01 . (8) Equation (8) can be interpreted as follows: The curva- ture of the optical surface is the difference between the projected curvatures of incoming and the outgoing wave fronts, which is a direct extension of Eq. (1). Projection essentially increases curvature in the sagittal direction and decreases it in the tangential direction. This phe- nomenon was outlined in an earlier paper.8 Given the curvature of the outgoing wave front, we can calculate the foci of the outgoing beam, g1and g2, and the irradiance produced at the target surface according to Eq. (4): Et IIn psps pg1 兩兩 pg2 tg1 兩兩 tg2 NtOut. (9) Here I(In) is the intensity of the source in the direction In, i.e., the irradiance produced by the source on the unit sphere in that direction. The ﬁrst fraction in Eq. (9) de- scribes the irradiance produced by the source at the posi- tion of the optical surface, measured in the tangent plane of the incoming radiation. The second fraction describes propagation of the outgoing beam characterized by the foci g1and g2between the optical surface at pand the target point t. Finally, the last scalar product accounts for the incidence angle with respect to the target surface normal Nt. Note that the positions of the foci g1and g2 refer to the outgoing wave front. Their relation to the curvature of the incoming radiation and that of the opti- cal surface is quantiﬁed by Eq. (5). If one chooses orthonormal coordinate systems aligned with a sagittal and a tangential plane on all three sur- faces, the projection operators are diagonal, which consid- erably simpliﬁes Eq. (5). But in this coordinate system, in general the curvature tensor of the optical surface is not diagonal. One might be tempted to restrict the cur- vature of the optical surface by requiring that the eigen- vectors of the curvature tensor also be aligned with the sagittal and tangential directions, deﬁned by the inci- dence direction. After all, the distance to these foci is all that matters for the irradiance. However, for reasons outlined below, this choice violates the integrability con- dition and therefore in general rules out a solution of the differential equations as a smooth, continuous optical sur- face. Equation (5) relates the curvature of the optical surface to its slope. The great achievement of Gauss was to prove that slope and curvature are inherent properties of a surface, which do not depend on a particular parameter- ization. We made use of this fact by choosing an abstract formulation in the previous derivations in order to better illustrate the underlying principles. If, however, one proceeds to actually construct an opti- cal surface, one cannot avoid resorting to a particular rep- resentation based on some coordinate system. As we have seen, it is not possible to diagonalize the projection and the curvature tensors simultaneously. Equations (1), (4), and (5) translate to a set of two elliptical nonlin- ear partial differential equations. Equivalent equations were derived by Monge and Ampe re more than 200 years ago.9The general problem for a point source was pre- sented in the form of an elliptical equation of the Monge Ampe re type by Komissarov and Boldyrev.10 The prob- lem that we have treated has also been analyzed in detail by Schruben.11 This author used an equivalent but nev- ertheless different notation and also showed that the problem is governed by an elliptic MongeAmpe re differ- ential equation. However, the earlier authors give no hint about actually generating a solution for a particular desired distribution. The reason for us to use a different formulation is that our primary goal is to actually ﬁnd an algorithm that generates a numerical solution in the gen- eral case. We did not attempt a closed-form solution be- cause we did not wish to restrict the design input to func- tions represented as closed-form expressions. Indeed, although the numerical problem is quite demanding and challenging, we were able to devise a numerical algorithm for ﬁnding a numerical solution following standard nu- merical procedures outlined in Ref. 12. On the basis of the examples we have tackled so far, we are conﬁdent that we have developed a powerful tool applicable to an ex- tremely large class of problems. Suppose we have integrated the second-order equation and determined the normal to the surface. Is it possible to determine the reﬂector surface? It has been shown13 that a necessary and sufﬁcient condition for the existence of a surface everywhere perpendicular to a vector ﬁeld N, i.e., for the existence of a quasipotential, is that this vec- tor ﬁeld satisfy NcurlN0. (10) This constitutes an additional condition to be imposed on the on the second derivative of the surface parameteriza- tion. An equivalent condition was found by Parkyn.14 One is free to choose boundary conditions for a particular solution; however, one cannot freely choose a complete boundary curve for the optical surface. In general, such a choice would ﬁx at least one component of the normal vector and therefore violate either Eq. (1) or, most prob- ably, the integrability condition. Therefore we chose to require the boundary of the optical surface to lie on an a priori deﬁned two-dimensional manifold such as a cyl- 592 J. Opt. Soc. Am. A /Vol. 19, No. 3 /March 2002 H. Ries and J. Muschaweck inder or a cone. We found that with this boundary con- dition the solution is unique up to a scaling factor. 3. NUMERICAL SOLUTION AND BOUNDARY CONDITIONS In order to determine the shape of the optical surface, we have to solve the set of differential equations. Let us re- view them in turn: Equation (1), which speciﬁes which point on the tar- get surface is illuminated by a given point on the optical surface. Equation (9), which connects the desired irradiance on the target to the focal points of the outgoing wave front. Equation (5), which connects the curvature of the outgoing wave front to the curvature of the incoming wave front and the curvature of the optical surface. Equation (10), the integrability condition. Equation (10) has a key role. It may seem that the de- sired irradiance distribution on the target may be achieved with an entire family of curvature tensors, and one may simply choose a convenient representative. In- deed, solutions based on such an arbitrary choice have been shown to work in particular, usually highly symmet- ric, conﬁgurations.15 However, we found that an arbi- trary choice in the general case violates the integrability condition. Therefore we have added the integrability condition explicitly to the set of determining equations. The set of equations is nonlinear and elliptical. As out- lined in Ref. 12, equations of this type need special care. We have been able to devise an algorithm by discretiza- tion of the optical surface to a grid of points and subse- quently identifying derivatives with ﬁnite differences. In addition to the differential equations, one can specify boundary conditions. The boundary conditions also need special consideration as indicated by Schruben.11 If we were simply to specify the outer boundary of the optical surface as a curve in space, this would also specify one tangent at the rim and therefore in general violate Eq. (1). Therefore we can require only that the boundary of the optical surface lie on some pre- deﬁned surface. This allows us, for example, to require that the boundary of the optics be on a cylinder but not that it be on a circle. In addition, we found that because in the limit of a small source the optical laws scale, we can require a cer- tain average size of the optical surface (or average dis- tance from the source). However, this condition may vio- late energy conservation, because the entire optical surface must intercept a total power from the source that precisely matches the integral of the desired irradiance on the target. An additional restriction applies to refractive optics. The maximum angle by which a ray can be redirected by a lens is limited by total internal reﬂection. When tailor- ing refractive optics, we did encounter occasionally un- physical solutions, which attempt to change the direction of the rays by too large an angle, meaning that in reality, total internal reﬂection would occur for part of or the en- tire optical surface. The two caustic surfaces are the reason for the exis- tence of three topologically distinct solutions. We have explained that the caustic surfaces must never intersect the target surface. Therefore in general there are dis- tinct solutions: one in which both caustic surfaces lie be- tween source and target, one in which both caustic sur- faces lie outside the region between source and target, and a third in which one caustic surface lies between source and target and one does not. Apart from these distinctions and the restrictions outlined above, we found that if a suitable set of boundary conditions is speciﬁed, a solution can be found by tailoring, and this solution ap- pears to be unique. Fig. 3. Sketch of the setting: The lens is made to collect light from an isotropic point source and redistribute it onto a target perpendicular to the axis. Length units are arbitrary; the draw- ing is not to scale. Fig. 4. Shape of the freeform lens that casts the letters OEC on a square background, as it would appear under oblique illumina- tion. A slight bulging is visible, which remotely resembles the letter shape. Note that the rim lies on a cone as required by the boundary condition. H. Ries and J. Muschaweck Vol. 19, No. 3/ March 2002 /J. Opt. Soc. Am. A 593 4. SOLVED EXAMPLES As a demanding example and proof of concept, we speci- ﬁed a desired irradiance distribution that is constant over a square target except for the contours of the letters ‘‘OEC,’’ for which we speciﬁed an intensity three times that of the background. All rays incident on the optical surface are to be redistributed on the target with no rays outside. For this case we designed a refractive (lens) sur- face to perform the task in a setting sketched in Fig. 3. Of course, it would be a simple task to create any de- sired irradiance distribution in this setting with a simple slide projector. However, such an approach wastes radia- tion by absorption in the slide. In contrast, the method that we outline here creates the desired irradiance distri- bution purely by redirecting the radiation. Thereby all intercepted radiation is redirected to the target. The shape of the lens is depicted in Figs. 4 and 5. By careful inspection, one has the impression of recognizing regions on the lens surface singled out by slightly differ- ent curvatures, which resemble the letters on the target although the shape is distorted. In retrospect, this is easy to understand since the regions of the lens contrib- uting to the letters have to be three times as large in area as a result of energy conservation. The rim of the lens was made to lie on a cone, oriented with its axis toward the center of the desired distribution and a tip half-angle of roughly 33 deg. The energy balance determines the exact value of the tip angle. We veriﬁed the optical performance of the tailored free- form lens by importing the shape into a commercial ray- tracing program16 and performed an independent evalua- tion by tracing 9 106rays from a point source through the lens and computing their intersection with the target plane. The result is shown in Fig. 6. Figure 7 plots the intensity along the horizontal cross section through the center in order to quantitatively verify the intensity ratio of zero (outside) to 1 to 3. 5. CONCLUSIONS We have presented a method for designing freeform opti- cal surfaces tailored to speciﬁc illumination tasks. Given a small source with known intensity distribution, given also a desired irradiance distribution on a target surface, we have devised a constructive procedure to design the shape of an optical surface, either reﬂective or refractive, that redistributes the radiation from the source onto the target, precisely matching the desired distribution. Tai- loring draws on principles of nonimaging optics, applied, as far as we know for the ﬁrst time, to problems without Fig. 5. Shape of the freeform lens in an isosag representation. Lines indicate equal sag. Spacing is 0.01 unit. The ragged rim is an artifact of the representation. The rim of the lens was made to lie on a cone oriented with its axis toward the center of the desired distribution. Fig. 6. Irradiance distribution produced by the freeform tailored lens based on Monte Carlo ray tracing of 9 106rays. Inter- section of the rays with a plane at the position of the nominal tar- get was binned by 129 129 pixels. Intensity variations are entirely compatible with shot noise. Fig. 7. Horizontal cross section through the irradiance distribu- tion in Fig. 5 at vertical position 0. Irradiance is in arbitrary units. Note the contrast of 0/1/3. Noise is completely within statistical bounds. 594 J. Opt. Soc. Am. A /Vol. 19, No. 3 /March 2002 H. Ries and J. Muschaweck any symmetry. Because of its generality this approach is very promising for a large variety of applications. We have presented idealized results in which we delib- erately excluded many real-world effects such as ﬁnite source size, Fresnel reﬂections and, above all, manufac- turing errors such as slope errors of the surface. We did this in order to show the strength and limitations of the mathematical procedure in the sense of a proof of concept. We found that the numerical procedure is extremely fast and efﬁcient. For the example presented, in which the surface was rendered as a grid of 129 129 points, the computation required roughly one hour on a standard contemporary personal computer. The numerical preci- sion was far beyond manufacturing capabilities. In fact, we found it more difﬁcult to assess the perfor- mance of the tailored optical surface with ray tracing than to design it. For an illustration, imagine an ambi- tious claim that within a target area of, say, 1 m2there is not a single squared millimeter on which the irradiance is more than 1% off the nominal value. In fact, we believe that apart from manufacturing tolerances, the pure com- putational procedure outlined here is quite capable of such outstanding performance because we have veriﬁed that the tolerance in selected outgoing rays is of the order of 1 m. To substantiate such a claim by a Monte Carlo ray-tracing technique, one would need more than 104rays per bin in order to reach a statistical error below 1%. On the other hand, to match the spatial resolution, one would need at least 106pixels. These requirements combine to a requirement of 1010 rays. This is beyond the capability of today’s computers, which typically require of the order of 103s per ray. We are well aware that our example is of no practical value. The reason we chose it was, again, to demonstrate a capability that other known design methods can hardly match. The desired irradiance distribution is what we would call rich in detail. It is hard to imagine that a pa- rameterization of a surface based on few parameters for such a task would be feasible. The example chosen deliberately involves a discontinu- ous desired irradiance distribution that features steps in irradiance by a factor of 3. Because irradiance on the target essentially is linked to the curvature of the tailored surface, the optical surface designed is continuous and smooth. Only the curvature exhibits steps. The tailored surface is geometrically smooth (G1) , even for disconti- nous desired intensity distribution. Since intensity is re- lated to the curvature of the optical surface, in general, if the desired irradiance is cncontinuous, the tailored opti- cal surface will feature Gn2geometric continuity. This facilitates manufacturing. ACKNOWLEDGMENTS We thank Wolfgang Spirkl and Andreas Timinger for valuable suggestions and discussions. We are grateful to the reviewer for pointing out Refs. 10, 11, and 13. Harald Ries can be reached by e-mail at harald.ries @oec.net. REFERENCES 1. J. M. Gordon and H. Ries, ‘‘Tailored edge-ray concentrators as ideal second stages for Fresnel reﬂectors,’’ Appl. Opt. 32, 22432251 (1993). 2. H. Ries and R. Winston, ‘‘Tailored edge-ray reﬂectors for il- lumination,’’ J. Opt. Soc. Am. A 11 , 12601264 (1994). 3. H. Ries and J. Gordon, ‘‘Double-tailored imaging concentra- tors,’’ in Nonimaging Optics: Maximum Efﬁciency Light Transfer V, R. Winston, ed., Proc. SPIE 3781, 129 134 (1999). 4. J. M. Gordon and H. Ries, ‘‘Tailord edge-ray concentrators (TERC’s): Approaching the thermodynamic limit with new ideal second stages for Fresnel reﬂectors,’’ presented at The Biennial Congress of the International Solar Energy Soci- ety, Budapest, August 1923, 1993. 5. R. P. Friedman, J. M. Gordon, and H. Ries, ‘‘Compact high- ﬂux two-stage solar collectors based on tailored edge-ray concentrators,’’ Sol. Energy 56,607615 (1996). 6. L. D. Landau and E. M. Lifschitz, Klassische Feldtheorie. Vol. 2 of Lehrbuch der Theoretischen Physik, G. Heber, Ed. (Akademie Verlag, Berlin, 1904). 7. H. Ries and W. Spirkl, ‘‘Caustic and its use in designing op- timal absorber shapes for 2D concentrators,’’ in Nonimag- ing Optics: Maximum Efﬁciency Light Transfer III,R. Winston, ed., Proc. SPIE 2538,29 (1995). 8. R. Zaibel, E. Dagan, J. Karni, and H. Ries, ‘‘An astigmatic corrected target-aligned heliostat for high concentration,’’ Sol. Energy Mater. Sol. Cells 37, 191202 (1995). 9. G. Monge, Application de l’analyse a la ge ´ome ´trie, 5th ed. Corrected and annotated by M. Liouville (Bachelier, Paris, 1850). 10. V. D. Komissarov and N. G. Boldyrev, ‘‘The foundations of calculating specular prismatic ﬁttings,’’ Trudy VEI 43,661 (1941). 11. J. S. Schruben, ‘‘Formulation of a reﬂector design problem for a lighting ﬁxture,’’ J. Opt. Soc. Am. 62, 14981501 (1972). 12. W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vet- terling, Numerical Recipes, 2nd ed. (Cambridge U. Press, New York, 1992). 13. P. Moon, D. E. Spencer, The Photic Field (MIT Press, Cam- bridge, Mass., 1981). 14. W. A. Parkyn, ‘‘The design of illumination lenses via extrin- sic differential geometry,’’ in International Optical Design Conference, K. P. Thompson and L. R. Gardner, eds., Proc. SPIE 3482,191193 (1998). 15. W. A. Parkyn, ‘‘Illuminating Lens designed via extrinsic dif- ferential geometry,’’ U.S. patent 5,924,788, July 20, 1999. 16. LightTools 3.1, 2001, Optical Research Associates, www.op- ticalres.com, 3280 East Foothill Boulevard, Pasadena, Cali- fornia 91107. H. Ries and J. Muschaweck Vol. 19, No. 3/ March 2002 /J. Opt. Soc. Am. A 595 ... In dO, the merit function can also be at the irradiance level, when optimizing in the image space rather in the intersection space, e.g. in Refs. [48], [13]. This requires integrating a sensor image I(p) from intersections p, given a pre-defined pixel size. ... Article Full-text available Computational imaging systems algorithmically post-process acquisition images either to reveal physical quantities of interest or to increase image quality, e.g., deblurring. Designing a computational imaging system requires co-design of optics and algorithms, and recently Deep Lens systems have been proposed in which both components are end-to-end designed using data-driven end-to-end training. However, progress on this exciting concept has so far been hampered by the lack of differentiable forward simulations for complex optical design spaces. Here, we introduce\mathbf{dO}$(DiffOptics) to provide derivative insights into the design pipeline to chain variable parameters and their gradients to an error metric through differential ray tracing. However, straightforward back-propagation of many millions of rays requires unaffordable device memory, and is not resolved by prior works.$\mathbf{dO}$alleviates this issue using two customized memory-efficient techniques: differentiable ray-surface intersection and adjoint back-propagation. Broad application examples demonstrate the versatility and flexibility of$\mathbf{dO}$, including classical lens designs in asphere, double-Gauss, and freeform, reverse engineering for metrology, and joint designs of optics-network in computational imaging applications. We believe$\mathbf{dO}\$ enables a radically new approach to computational imaging system designs and relevant research domains.
... Non-imaging free-form lenses are formulated as an inverse problem with the aim to redistribute a light source with a known intensity profile into a desired output intensity [21]. A variety of different strategies have been implemented to tackle this inverse problem [21,[25][26][27][28][29], but a particularly elegant and general way to do it is through an optimal transport formulation [30,31]. This method makes no a priori symmetry assumptions on the lens shape and has been demonstrated to work for a variety of incident light distributions and output illumination patterns [32]. ...
Preprint
Full-text available
Metasurfaces provide a versatile platform for realizing ultrathin flat optics for use in a wide variety of optical applications. The design process involves defining or calculating the phase profile of the metasurface that will yield the desired optical output. Here, we present an inverse design method for determining the phase profile for shaping the intensity profile of a collimated incident beam. The model is based on the concept of optimal transport from non-imaging optics and enables a collimated beam with an arbitrary intensity profile to be redistributed to a desired output intensity profile. We derive the model from the generalized law of refraction and numerically solve the resulting differential equation using a finite-difference scheme. Through a variety of examples, we show that our approach accommodates a range of different input and output intensity profiles, and discuss its feasibility as a design platform for non-imaging optics.
... Tailoring design method is also used to develop illumination system design, and to build 3D freeform surfaces by solving nonlinear PDEs. The relationship between light vectors such as incident, normal, and exit, as well as constraints among the source of light with the irradiance distribution at the reference surface, are taken into account by the formed freeform optical surface [79]. ...
Article
Full-text available
Freeform optics has become the most prominent element of the optics industry. Advanced freeform optical designs supplementary to ultra-precision manufacturing and metrology techniques have upgraded the lifestyle, thinking, and observing power of existing humans. Imaginations related to space explorations, portability, accessibility have also witnessed sensible in today's time with freeform optics. Present-day design methods and fabrications techniques applicable in the development of freeform optics and the market requirements are focussed and explained with the help of traditional and non-traditional optical applications. Over the years, significant research is performed in the emerging field of freeform optics, but no standards are established yet in terms of tolerances and definitions. We critically review the optical design methods for freeform optics considering the image forming and non-image forming applications. Numerous subtractive manufacturing technologies including figure correction methods and metrology have been developed to fabricate extreme modern freeform optics to satisfy the demands of various applications such as space, astronomy, earth science, defence, biomedical, material processing, surveillance, and many more. We described a variety of advanced technologies in manufacturing and metrology for novel freeform optics. Next, we also covered the manufacturing-oriented design scheme for advanced optics. We conclude this review with an outlook on the future of freeform optics design, manufacturing and metrology.
... The design process often involves cumbersome calculations in order to determine the surface needed to project the desired intensity profile. Indeed, these geometries are typically defined as sum of spherical or aspherical lenses [151], through Q-polynomials description [152] or nonlinear partial differential equations [153]. Their fabrication is also a demanding task, in fact, manufacturing processes with high degrees of freedom are required [154], such as grinding, polishing [155], and ultra-precision turning [156]. ...
Thesis
Full-text available
In this thesis I will explain the properties of intersubband polaritons by a theoretical and experimental point of view. After a brief introduction to explain some general aspects of the intersubband polariton generation I will describe an experimental campaign with the aim to verify the possibility to create a polariton laser. The emission of light from a polaritonic state lacks a clear experimental verification, therefore this experiment may be the first observation of a polariton laser. The possibility to use graphene to design a better optical resonance will be largely discussed. In the end I will illustrate a concept to manipulate the light with 3D printed windows. The final objective being the design of windows able to remove aberration from laser source (in particular TeraHertz sources).
Article
Wide-field-of-view (FoV) Offner imaging spectrometers with freeform surfaces have been studied extensively in recent years. However, a design result with a large numerical aperture (NA) cannot be simultaneously obtained with this layout. We present the concept of a limited system in the tangential direction. Based on this insight, we present a new design method, to the best of our knowledge, based on the decenter anamorphic stop, which can achieve large NA in compact wide-FoV Offner imaging spectrometers with freeform surfaces. Compared to conventional imaging spectrometers with the same parameters, the light-gathering capacity of the decenter anamorphic stop-based imaging spectrometer is increased by more than 40%. In addition, based on the presented method, we design a compact imaging spectrometer with a wide FoV and large NA. The designed imaging spectrometer with a freeform surface has excellent performance. Finally, we fabricate and measure the freeform mirror. The surface irregularity of the freeform mirror is better than 1 / 30 λ ( λ = 632.8 n m ). The result shows that the Offner imaging spectrometer with a freeform surface can be fabricated and will play a significant role in the fields of aeronautical and astronautical remote sensing.
Article
Freeform tailoring optics are core components in advanced illumination and laser technologies. In this work, we propose an efficient and streamlined optics design method to face the challenges in its implementation. We integrate ray-mapping least squares calculations with freeform surface optimization tools embedded in commercial optical design software. As a result, the complexity of freeform illumination and laser beam shaping design is significantly reduced by only adapting a built-in merit function. We demonstrate the feasibility of our approach through three design examples that provide tailored irradiance distributions. Furthermore, the results highlight the efficiency, simplicity and versatility of our approach that can be readily applied to multiple configurations.
Article
A key challenge in tailoring compact and high-performance illumination lenses for extended non-Lambertian sources is to take both the étendue and the radiance distribution of an extended non-Lambertian source into account when redirecting the light rays from the source. We develop a direct method to tailor high-performance illumination lenses with prescribed irradiance properties for extended non-Lambertian sources. A relationship between the irradiance distribution on a given observation plane and the radiance distribution of the non-Lambertian source is established. Both edge rays and internal rays emanating from the extended light source are considered in the numerical calculation of lens profiles. Three examples are given to illustrate the effectiveness and characteristics of the proposed method. The results show that the proposed method can yield compact and high-performance illumination systems in both the near field and far field.
Article
A freeform reflector for illumination in a space-based remote sensor calibration system is designed to produce the desired patterns and improve efficiency. Owing to the excellent cosine properties and light source stability of the integrating sphere, it is widely used as the calibration source in optical remote sensing instruments. However, the light source cannot obtain uniform illumination on a large scale. Therefore, we propose a freeform reflector design for an integral sphere source, which generates an independent irradiance distribution and achieves large-scale illumination by arraying. The irradiance uniformity and efficiency obtained were greater than 98% and 90%, respectively.
Preprint
Full-text available
Adoption of LED lighting products have many advantages such as energy conservation and environmental protection. But the luminescent characteristics of LEDs cannot provide satisfactory lighting solutions to meet the requirements of today’s sophisticated applications. Freeform optics constitutes a new technique for regulating the energy output of LEDs into a customized form. However, many freeform optical elements are designed to be bulky since the point-like sources assumption. The designs that consider the spatial extents of LEDs usually have poorly controlled light patterns or low efficiencies. Here, we propose a very robust design method based on spatially variant convolution for 4D light field tailoring, which enables high-performing freeform illumination optics with small sizes. The proposed method is holistic that is suitable for various kinds of design requirements with no limitations of far-field approximation and Lambertian luminescent property. Several design examples including fabricated lenses and experimental tests are presented. The results realize compact, low-cost, and highly efficient freeform optics with customized lighting solutions. We believe that the proposed work can promote the research on enhancing the emission quality of LEDs and generating artificial light for application-specific lighting systems.
Article
Nonimaging optical design aiming at energy control has a wide range of applications in optoelectronics. A nonimaging optical system is composed of a light source, optical components, and a target screen, and can be described by an equation named light taming equation(LTE). Given the light source and prescribed target spot, the required freeform surfaces of the optical component can be obtained by solving the LTE. If the light source profile does not change, the optical surface has to make some suitable morphs when the target spot translates on the screen, and the group theory can describe these morph operators well. The basic LTE was established for a normal nonimaging optical system, which is to design an optical element for redirecting the lights from the source so that a prescribed light distribution is generated on a given target. A translation light taming equation(T-LTE) was derived for the case of only spot translating on the target screen, and an optical translation group(OTG) was introduced for describing all of the morph operators of the optical surface caused by light spot translation. There are multiple solutions for the same T-LTE, but the uniqueness of the T-LTE solution is necessary for OTG. Fortunately, the eikonal-energy(KE) mapping method can guarantee the uniqueness of the T-LTE solution, where K is the optical path length. The supporting quadric method(SQM) is one of the KE mapping methods when the nonimaging optical system has only one optical surface to be resolved. The LTE with SQM was deduced, and the OTG can be discussed in K-space. A deep neural network(DNN) was introduced to fit the KE mapping and spot translating operators to obtain the required optical surface. Taking the uniform square spot as an example, the SQM generated the sample data of spot translation to train DNN. The optical simulation results show that the error between the light distribution generated by the DNN and the standard uniform square spot is small, all within the order of 10⁻³, which indicates that the DNN and KE mapping method successfully realized the function of the OTG. There is a guiding significance for intelligent nonimaging optical design.
Article
Full-text available
Conventional heliostats suffer from astigmatism for non-normal incidence. For tangential rays the focal length is shortened while for sagittal rays it is longer than the nominal focal length. Due to this astigmatism it is impossible to produce a sharp image of the sun, and the rays will be spread over a larger area. In order to correct this the heliostat should have different curvature radii along the sagittal and tangential direction in the heliostat plane just like a non axial part of a paraboloid. In conventional heliostats, where the first axis, fixed with respect to the ground, is vertical while the second, fixed with respect to the reflector surface, is horizontal such an astigmatism correction is not practical because the sagittal and tangential directions rotate with respect to the reflector. We suggest an alternative mount where the first axis is oriented towards the target. The second axis, perpendicular to the first and tangent to the reflector, coincides with the tangential direction. With this mounting sagittal and tangential direction are fixed with respect to the reflector during operation. Therefore a partial astigmatism compensation is possible. We calculate the optimum correction and show the performance of the heliostat. We also show predicted yearly average concentrations.
Article
Full-text available
The edge-ray principle can be used to tailor a reflector. However, one set of edge rays already fully determines the reflector profile. We present a design method for tailoring compact compound elliptical concentrator (CEC)-type reflectors to a given source and a desired angular power distribution. Two reflected images of the source, one on each side of the source, contribute together with the direct radiation from the source to produce the desired power distribution. We determine the reflector profile by numerically solving a differential equation. No optimization is required. Beyond the angular region in which the power distribution can be strictly controlled, the power drops to zero in a finite decay range. This decay range becomes narrower as the reflector increases in size. We show a reflector for producing a strictly constant irradiance from −43 to 43 deg from a cylindrical source of constant brightness. The reflector extends to a maximum distance of 8 source diameters. No power is radiated beyond ± 50 deg.
Article
The caustic of a set of edge rays is defined as the set of intersection points of adjacent edge rays. For a body having a smooth differentiable contour, the caustic of its edge rays coincides with the contour of the body. Therefore one would assume that by calculating the caustic of the edge rays as they are produced by a 2D concentrator such as a trough, the optimal shape for the absorber, e.g. the minimal surface absorber capable of intercepting all rays, should also coincide with the shape of the caustic. We show that this conjecture is not valid in general, but only if the caustic indeed forms a closed smooth curve. For parabolic trough systems, the caustic intersects and forms closed domains for half rim angles of around 60 degrees and 120 degrees. In both cases the contour is not smooth. Therefore the optimal shape is not given by the domain enclosed by the caustic. We present a general recipe of how to construct minimum surface absorbers for given caustics in 2D and apply this to the case of trough parabolic concentrators. We show practical absorber shapes for parabolic troughs with various rim angles. The optimal contour depends discontinuously on the rim angle. The area of the optimum shape for a rim angle of 90 degrees is 0.72 of the area of the smallest cylindric absorber capable of intersecting all rays.
Article
The differential geometry of surfaces describes both intrinsic and extrinsic properties: the intrinsic ones, such as the metric, are inherent, whereas extrinsic ones describe shape relative to 3D space, such as the distribution of surface normal vectors. This paper describes how the surface normals of a non-rotationally-symmetric illumination lens are generated, via Snell's law in vector form, from (1) the distribution of light from a sources, and (2) a desired directional distribution of light exiting the lens. These two distributions are expressed as grids with cells of variable size but constant photometric flux, on the Gaussian sphere of directions. The grids must be sufficiently fine so as to generate enough surface normal vectors for accurate numerical generation of the requisite surface. The grids must have the same number of cells, and the same topology (i.e., rectangular vs. polar). The source-intensity grid must be adjusted to account for Fresnel reflection losses. For an unfaceted (smooth-surfaced) lens, the array of surface normal vectors must be adjusted for equality of the crossed partial derivatives. This class of lenses has only recently become producible due to the advent of electric- discharge machining for the shaping of non-rotationally symmetric injection molds for plastic lenses.
Article
A novel method is described of designing lenses for general illumination tasks that are not circularly symmetric. The crux of the method is the specification of the illumination task by a grid on the unit sphere of directions. This grid, or tessellation, has cells which vary in solid angle such that each encompasses the same luminous flux: high intensity corresponds to small cells, and vice versa. Another grid, having the same topology and number of cells, is formed according to the intensity distribution of the source. An illumination lens must then transform the source distribution into the task distribution, via one or more refractions. Thus the direction vector of a cell in the source grid must be redirected into that of the corresponding cell in the task grid. Snell's law in vector form enables the derivation of a corresponding surface normal vector, or sequence of normal vectors, that will accomplish this redirection. Extrinsic differential geometry is then used to generate a lens surface having, as closely as possible, this distribution of surface normals, which must be irrotational to generate a smooth surface. This class of lenses has only recently become producible due to the advent of electric-discharge machining for the shaping of non-rotationally symmetric injection molds for plastic lenses.
Article
We present a new approach in optical design whereby two- stage axisymmetric reflectors are tailored with a completely imaging strategy, and can closely approach the thermodynamic limit to radiation concentration at near-maximum collection efficiency. Practical virtues include: (1) an inherent large gap between the receiver and the second-stage mirror; (2) an upward-facing receiver; (3) the possibility of compact units (large rim angles), i.e., low ratios of total depth to total width; and (4) no chromatic aberration. We describe how one can tailor both the primary and secondary mirrors so as to insure that spherical aberration is eliminated in all orders, and circular coma is canceled up to first order in the angle subtended by the radiation source. An illustrative solution that attains about 93% of the thermodynamic limit to concentration is presented for a far-field source, as is common in solar energy and infrared detection applications. Double-tailored imaging concentrators are similar in principle to complementary Cassegrain concentrators that comprise a paraboloidal primary mirror and the inner concave surface of a hyperboloid secondary reflector, but have monotonic contours that are substantially different with far superior flux concentration.
Article
The problem of designing a reflector to distribute the illumination of a nonisotropic point source on a plane aperture according to a pre-assigned pattern is analyzed. An integral equation and equivalent partial differential equation are derived. The form of the latter reveals this reflector-design problem to be a singular elliptic Monge—Ampère boundary-value problem.


text:
# All Questions

29 views

### Can a flow of current be produced in a wire if its is bombarded with Electromagnetic waves?

If you have a conducting wire and you bombarded it with EM waves, is there a frequency (lower then the threshold frequency to liberate the electrons) that would induce a current within the wire. I ...
9 views

### Assumptions that led to Wigner surmise of nuclear level spacings

In 1957, Wigner surmised (guessed) that the probability density of the spacing between adjacent energy levels of complex nuclei is given by $$P(s) = \frac{\pi}{2}s e^{\frac{-\pi s^2}{4}}$$ This is ...
28 views

### The question about an observer and an observed environment [on hold]

Imagine that Bob is observing everything that is not him (environment). The experiment: the environment changes somewhere far away from Bob. The conclusion: The change propagates in the ...
12 views

### the total deviation produced by slab when Ray comes out [on hold]

Q) A ray of light travelling in air is incident at grazing incidence on a slab with refractive index,$n(y)=[ky^{3/2}+1]^{1/2}$ where $k= 1m^{-3/2}$ and follows path as shown in figure. What is the ...
27 views

### Homework question on the Poynting vector

I need to show what the Poynting vector becomes when i insert the definition of A through creation and annihilation operators, what i get is the expected answer plus a term containing two creation ...
21 views

### How many times does baryogenesis happen during the history of the universe?

Does the baryogenesis happen both in GUT epoch and electroweak Epoch during the history of the universe? If yes, are there any differences during this two times?
48 views

### What is the physical implication of a homomorphism between SU(2) and SO(3)?

It can be shown that there exists a homomorphism between $SU(2)$ and $SO(3)$. What is the physical implication of a homomorphism actually? I know that in physics $SU(2)$ acts on a space of spin 1/2 ...
35 views

### Metal density increase [on hold]

Can you increase the density of a material (metal) by surface treatment? For example commercial copper. If we take the ordinary OF copper (99.95%), I would expect not to have the tabulated 8.92/8.96 ...
21 views

### Background and references for topological string theory

Where would one start to learn about topological string theory knowing no string theory? What are the general prerequisites and what kind of background in QFT is needed? Is a conventional book like ...
25 views

### What causes a charged particle to follow circular motion in a magnetic field? [duplicate]

I would like to know why particles follow a circular path in a magnetic field, what path they trace and how this path is measured and controlled?
25 views

### Determine the shear in dislocation [on hold]

I don't know how to solve this kind of problem: Given a crystal and its volume V. This Crystal is subjected to a shear stress making a slip plane (h k l ) in a slip direction [h k l]. Given also the ...
68 views

### How do light waves get their size?

An atom or (small) molecule has the size of about 100pm. Elektromagnetic waves range from about 0.1nm up to 1 km. The most common way waves (like light) are caused by 'jumping' electrons to another ...
4k views

### Supergun Launching of Satellites

I should say first that I don't believe this is a feasible launch method, otherwise NASA and other space agencies would be using it by now. It's based on this BBC news story Saddam Hussein's Supergun ...
20 views

### The direction of frictional force in circular turning

Why does the frictional force in case of circular motion point towards the center even though the motion is tangential to the radius?
14 views

### Relationship between refractive index and concentration [on hold]

what does the gradient of the graph of refractive index against solution concentration mean
63 views

### Conformal Gravity

Lubos, in his comment to a question, says that (http://physics.stackexchange.com/q/61281) First of all, one can't gauge a symmetry without modifying (enriching) the field contents. Gauging a ...
18 views

### What is the name of basis states of bulk k.p Hamiltonian?

A k.p Hamiltonian for a bulk material can be represented by 8x8 matrix in basis of $|S\uparrow\rangle$, $|S\downarrow\rangle$, $|X\uparrow\rangle$, $|Y\uparrow\rangle$, $|Z\uparrow\rangle$, ...
15 views

### What is the vapour pressure of carbon at room temperature?

The vapour pressure of carbon at 3642°C is one atmosphere (this is the sublimation temperature given in Wikipedia). At lower temperatures, the vapour pressure will be less, but of course it will ...
44 views

### Why $U(1)_Y$ hypercharge rather than $U(1)_\text{em}$ electromagnetism?

In the Standard Model we have $SU(2)_I\times U(1)_Y$, where $U(1)_Y$ is weak hypercharge and $SU(2)_I$ is the symmetry group of weak isospin. Why do we introduce $U(1)_Y$ of weak hypercharge rather ...
26 views

### Working out total electrostatic force [on hold]

To work out this question, would you find the electrostatic forces between z and y and z and x using the formula F=kq1q2/d^2, and just add them?
26 views

### Can vernier scale ratio be less than 1?

I came across a question (See Image) that states about two vernier scales The first vernier scale is the top-most figure in the image: It has 10 vernier scale division equal to 9 main scale division ...
79 views

35 views

### Barycenter correction for inner planets in solar system?

I want to run simulations of the solar system, and I've read in some papers that to speed up the simulations it is useful to use the barycenter correction of the inner planets (from Mercury to Mars), ...
27 views

### BTZ Black holes

I computed the two point function for two scalar fields in BTZ black hole that is defined as a local $AdS_3$ space time with discrete identification, as defined in many papers. Referring to this ...
44 views

### Can air be considered incompressible as long flow velocities are less than 100 m/s (Ma = 0.3)?

Multiple sources state that steady air flow (in open systems) can be considered incompressible at velocities less than 100 m/s (Ma = 0.3). Example: Deborah A. Kaminski, Michael K. Jensen, ...
7k views

### Why would a spinning space station create a centrifugal force on an astronaut rather than simply spinning around him/her?

We often see films with spinning space station that create artificial gravity by having the astronauts pulled outwards by centrifugal force. I'd like to know if this would really happen, and if so, ...
70 views

38 views

### Global and relative phases of kets in QM

In one of the questions I'm trying to solve it is asked to, first, compute probabilities for the respective results of the Stern-Gerlach measurements performed on each state $\lvert\psi_1\rangle$, ...
815 views

### Orbits around the Photon sphere of a black hole (Schwarzschild coordinates)

This is a follow-up question to the answer given at What is the exact gravitational force between two masses including relativistic effects?. Unfortunately the author hasn't been online for a few ...
28 views

### Temperature of a single dipole [duplicate]

When taught about the temperature of a system from an illustrative point of view, I was told to think of molecules vibrating and hitting one another. If the molecules vibrate faster and faster (and in ...
82 views
+100

### Thermodynamic equilibrium - will the piston move?

Friend asked me this question and I didn't manage to solve it with basic thermodynamic reasoning. I believe this problem is easly solvable via numeric methods by choosing specific systems, though I ...

format: plain


text:
# Uniform distribution (continuous)

In probability theory and statistics, the continuous uniform distribution or rectangular distribution is a family of symmetric probability distributions. The distribution describes an experiment where there is an arbitrary outcome that lies between certain bounds.[1] The bounds are defined by the parameters, a and b, which are the minimum and maximum values. The interval can be either closed (eg. [a, b]) or open (eg. (a, b)).[2] Therefore, the distribution is often abbreviated U (a, b), where U stands for uniform distribution.[1] The difference between the bounds defines the interval length; all intervals of the same length on the distribution's support are equally probable. It is the maximum entropy probability distribution for a random variable X under no constraint other than that it is contained in the distribution's support.[3]

Notation Probability density functionUsing maximum convention Cumulative distribution function ${\displaystyle {\mathcal {U}}(a,b)}$ or ${\displaystyle \mathrm {unif} (a,b)}$ ${\displaystyle -\infty ${\displaystyle x\in [a,b]}$ ${\displaystyle {\begin{cases}{\frac {1}{b-a}}&{\text{for }}x\in [a,b]\\0&{\text{otherwise}}\end{cases}}}$ ${\displaystyle {\begin{cases}0&{\text{for }}x ${\displaystyle {\tfrac {1}{2}}(a+b)}$ ${\displaystyle {\tfrac {1}{2}}(a+b)}$ any value in ${\displaystyle (a,b)}$ ${\displaystyle {\tfrac {1}{12}}(b-a)^{2}}$ 0 ${\displaystyle -{\tfrac {6}{5}}}$ ${\displaystyle \ln(b-a)\,}$ ${\displaystyle {\begin{cases}{\frac {\mathrm {e} ^{tb}-\mathrm {e} ^{ta}}{t(b-a)}}&{\text{for }}t\neq 0\\1&{\text{for }}t=0\end{cases}}}$ ${\displaystyle {\begin{cases}{\frac {\mathrm {e} ^{itb}-\mathrm {e} ^{ita}}{it(b-a)}}&{\text{for }}t\neq 0\\1&{\text{for }}t=0\end{cases}}}$

## Definitions

### Probability density function

The probability density function of the continuous uniform distribution is:

${\displaystyle f(x)={\begin{cases}{\frac {1}{b-a}}&\mathrm {for} \ a\leq x\leq b,\\[8pt]0&\mathrm {for} \ xb\end{cases}}}$

The values of f(x) at the two boundaries a and b are usually unimportant because they do not alter the values of the integrals of f(x) dx over any interval, nor of x f(x) dx or any higher moment. Sometimes they are chosen to be zero, and sometimes chosen to be 1/(b  a). The latter is appropriate in the context of estimation by the method of maximum likelihood. In the context of Fourier analysis, one may take the value of f(a) or f(b) to be 1/(2(b  a)), since then the inverse transform of many integral transforms of this uniform function will yield back the function itself, rather than a function which is equal "almost everywhere", i.e. except on a set of points with zero measure. Also, it is consistent with the sign function which has no such ambiguity.

Graphically, the probability density function is portrayed as a rectangle where (b-a) is the base of the rectangle and (1/(b-a)) is the height of the rectangle. As the distance between a and b increases, the density at any particular value within the distribution boundaries decreases.[4] Since the probability density function integrates to 1, the height of the probability density function decreases as the base length increases.[4]

In terms of mean μ and variance σ2, the probability density may be written as:

${\displaystyle f(x)={\begin{cases}{\frac {1}{2\sigma {\sqrt {3}}}}&{\mbox{for }}-\sigma {\sqrt {3}}\leq x-\mu \leq \sigma {\sqrt {3}}\\0&{\text{otherwise}}\end{cases}}}$

#### Example 1. Using the Uniform Probability Density Function[5]

For random variable X

X~U(0,23)

Find P(2 < X < 18):

P(2 < X < 18) = (18-2)*(1/(23-0)) = 16/23.

In graphical representation of uniform distribution function [f(x) vs x], the area under the curve within the specified bounds displays the probability (shaded area is depicted as a rectangle). For this specific example above, the base would be (18-2) and the height would be (1/23).[5]

#### Example 2. Using the Uniform Probability Density Function (Conditional)[5]

For random variable X

X~U(0,23)

Find P(12 < X | X > 18):

P(X > 12 | X > 8) = (23-12)*(1/(23-8))=11/15.

The example above is for a conditional probability case for the uniform distribution: given X > 8 is true, what is the probability that X > 12. Conditional probability changes the sample space so a new interval length (b-a) has to be calculated, where b is 23 and a is 8.[5] The graphical representation would still follow Example 1, where the area under the curve within the specified bounds displays the probability and the base of the rectangle would be (23-12) and the height (1/15).[5]

### Cumulative distribution function

${\displaystyle F(x)={\begin{cases}0&{\text{for }}xb\end{cases}}}$

Its inverse is:

${\displaystyle F^{-1}(p)=a+p(b-a)\,\,{\text{ for }}0

In mean and variance notation, the cumulative distribution function is:

${\displaystyle F(x)={\begin{cases}0&{\text{for }}x-\mu <-\sigma {\sqrt {3}}\\{\frac {1}{2}}\left({\frac {x-\mu }{\sigma {\sqrt {3}}}}+1\right)&{\text{for }}-\sigma {\sqrt {3}}\leq x-\mu <\sigma {\sqrt {3}}\\1&{\text{for }}x-\mu \geq \sigma {\sqrt {3}}\end{cases}}}$

and the inverse is:

${\displaystyle F^{-1}(p)=\sigma {\sqrt {3}}(2p-1)+\mu \,\,{\text{ for }}0\leq p\leq 1}$

### Generating functions

#### Moment-generating function

${\displaystyle M_{x}=E(e^{tx})={\frac {e^{tb}-e^{ta}}{t(b-a)}}\,\!}$

from which we may calculate the raw moments m k

${\displaystyle m_{1}={\frac {a+b}{2}},\,\!}$
${\displaystyle m_{2}={\frac {a^{2}+ab+b^{2}}{3}},\,\!}$
${\displaystyle m_{k}={\frac {1}{k+1}}\sum _{i=0}^{k}a^{i}b^{k-i}.\,\!}$

For the special case a = b, that is, for

${\displaystyle f(x)={\begin{cases}{\frac {1}{2b}}&{\text{for}}\ -b\leq x\leq b,\\[8pt]0&{\text{otherwise}},\end{cases}}}$

the moment-generating functions reduces to the simple form

${\displaystyle M_{x}={\frac {\sinh bt}{bt}}.}$

For a random variable following this distribution, the expected value is then m1 = (a + b)/2 and the variance is m2  m12 = (b  a)2/12.

#### Cumulant-generating function

For n  2, the nth cumulant of the uniform distribution on the interval [-1/2, 1/2] is Bn/n, where Bn is the nth Bernoulli number.[7]

### Standard uniform

Restricting ${\displaystyle a=0}$ and ${\displaystyle b=1}$, the resulting distribution U(0,1) is called a standard uniform distribution.

One interesting property of the standard uniform distribution is that if u1 has a standard uniform distribution, then so does 1-u1. This property can be used for generating antithetic variates, among other things. In other words, this property is known as the inversion method where the continuous standard uniform distribution can be used to generate random numbers for any other continuous distribution.[4] If u is a uniform random number with standard uniform distribution (0,1), then x = Inverse of F(u) generates a random number x from any continuous distribution with the specified cumulative distribution function F.[4]

### Relationship to other functions

As long as the same conventions are followed at the transition points, the probability density function may also be expressed in terms of the Heaviside step function:

${\displaystyle f(x)={\frac {\operatorname {H} (x-a)-\operatorname {H} (x-b)}{b-a}},\,\!}$

or in terms of the rectangle function

${\displaystyle f(x)={\frac {1}{b-a}}\,\operatorname {rect} \left({\frac {x-\left({\frac {a+b}{2}}\right)}{b-a}}\right).}$

There is no ambiguity at the transition point of the sign function. Using the half-maximum convention at the transition points, the uniform distribution may be expressed in terms of the sign function as:

${\displaystyle f(x)={\frac {\operatorname {sgn} {(x-a)}-\operatorname {sgn} {(x-b)}}{2(b-a)}}.}$

## Properties

### Moments

The mean (first moment) of the distribution is:

${\displaystyle E(X)={\frac {1}{2}}(a+b).}$

The second moment of the distribution is:

${\displaystyle E(X^{2})={\frac {1}{3}}(a^{2}+ab+b^{2}).}$

In general, the n-th moment of the uniform distribution is:

${\displaystyle E(X^{n})={\frac {1}{n+1}}\sum _{k=0}^{n}a^{k}b^{n-k}.}$

The variance (second central moment) is:

${\displaystyle V(X)={\frac {1}{12}}(b-a)^{2}}$

### Order statistics

Let X1, ..., Xn be an i.i.d. sample from U(0,1). Let X(k) be the kth order statistic from this sample. Then the probability distribution of X(k) is a Beta distribution with parameters k and nk + 1. The expected value is

${\displaystyle \operatorname {E} (X_{(k)})={k \over n+1}.}$

This fact is useful when making Q–Q plots.

The variances are

${\displaystyle \operatorname {V} (X_{(k)})={k(n-k+1) \over (n+1)^{2}(n+2)}.}$

### Uniformity

The probability that a uniformly distributed random variable falls within any interval of fixed length is independent of the location of the interval itself (but it is dependent on the interval size), so long as the interval is contained in the distribution's support.

To see this, if X ~ U(a,b) and [x, x+d] is a subinterval of [a,b] with fixed d > 0, then

${\displaystyle P\left(X\in \left[x,x+d\right]\right)=\int _{x}^{x+d}{\frac {\mathrm {d} y}{b-a}}\,={\frac {d}{b-a}}\,\!}$ which is independent of x. This fact motivates the distribution's name.

### Generalization to Borel sets

This distribution can be generalized to more complicated sets than intervals. If S is a Borel set of positive, finite measure, the uniform probability distribution on S can be specified by defining the pdf to be zero outside S and constantly equal to 1/K on S, where K is the Lebesgue measure of S.

## Statistical Inference

### Estimation of parameters

#### Estimation of maximum

##### Minimum-variance unbiased estimator

Given a uniform distribution on [0, b] with unknown b, the minimum-variance unbiased estimator (UMVUE) for the maximum is given by

${\displaystyle {\hat {b}}_{\text{UMVU}}={\frac {k+1}{k}}m=m+{\frac {m}{k}}}$

where m is the sample maximum and k is the sample size, sampling without replacement (though this distinction almost surely makes no difference for a continuous distribution). This follows for the same reasons as estimation for the discrete distribution, and can be seen as a very simple case of maximum spacing estimation. This problem is commonly known as the German tank problem, due to application of maximum estimation to estimates of German tank production during World War II.

##### Maximum likelihood estimator

The maximum likelihood estimator is given by:

${\displaystyle {\hat {b}}_{ML}=m}$

where m is the sample maximum, also denoted as ${\displaystyle m=X_{(n)}}$ the maximum order statistic of the sample.

##### Method of moment estimator

The method of moments estimator is given by:

${\displaystyle {\hat {b}}_{MM}=2{\bar {X}}}$

where ${\displaystyle {\bar {X}}}$ is the sample mean.

#### Estimation of midpoint

The midpoint of the distribution (a + b) / 2 is both the mean and the median of the uniform distribution. Although both the sample mean and the sample median are unbiased estimators of the midpoint, neither is as efficient as the sample mid-range, i.e. the arithmetic mean of the sample maximum and the sample minimum, which is the UMVU estimator of the midpoint (and also the maximum likelihood estimate).

### Confidence interval

#### For the maximum

Let X1, X2, X3, ..., Xn be a sample from U( 0, L ) where L is the population maximum. Then X(n) = max( X1, X2, X3, ..., Xn ) has the density[8]

${\displaystyle f_{n}(X_{(n)})=n{\frac {1}{L}}\left({\frac {X_{(n)}}{L}}\right)^{n-1}=n{\frac {X_{(n)}^{n-1}}{L^{n}}},0

The confidence interval for the estimated population maximum is then ( X(n), X(n) / α1/n ) where 100(1  α)% is the confidence level sought. In symbols

${\displaystyle X_{(n)}\leq L\leq X_{(n)}/\alpha ^{1/n}}$

### Hypothesis Testing

In statistics, when a p-value is used as a test statistic for a simple null hypothesis, and the distribution of the test statistic is continuous, then the p-value is uniformly distributed between 0 and 1 if the null hypothesis is true.

## Occurrence and applications

The probabilities for uniform distribution function are simple to calculate due to the simplicity of the function form.[2] Therefore, there are various applications that this distribution can be used for as shown below: hypothesis testing situations, random sampling cases, finance, etc. Furthermore, generally, experiments of physical origin follow a uniform distribution (eg. emission of radioactive particles).[1] However, it is important to note that in any application, there is the unchanging assumption that the probability of falling in an interval of fixed length is constant.[2]

### Economics Example for uniform distribution

In the field of economics, usually demand and replenishment may not follow the expected normal distribution. As a result, other distribution models are used to better predict probabilities and trends such as Bernoulli process.[9] But according to Wanke (2008), in the particular case of investigating lead-time for inventory management at the beginning of the life cycle when a completely new product is being analyzed, the uniform distribution proves to be more useful.[9] In this situation, other distribution may not be viable since there is no existing data on the new product or that the demand history is unavailable so there isn't really an appropriate or known distribution.[9] The uniform distribution would be ideal in this situation since the random variable of lead-time (related to demand) is unknown for the new product but the results are likely to range between a plausible range of two values.[9] The lead-time would thus represent the random variable. From the uniform distribution model, other factors related to lead-time were able to be calculated such as cycle service level and shortage per cycle. It was also noted that the uniform distribution was also used due to the simplicity of the calculations.[9]

### Sampling from an arbitrary distribution

The uniform distribution is useful for sampling from arbitrary distributions. A general method is the inverse transform sampling method, which uses the cumulative distribution function (CDF) of the target random variable. This method is very useful in theoretical work. Since simulations using this method require inverting the CDF of the target variable, alternative methods have been devised for the cases where the cdf is not known in closed form. One such method is rejection sampling.

The normal distribution is an important example where the inverse transform method is not efficient. However, there is an exact method, the Box–Muller transformation, which uses the inverse transform to convert two independent uniform random variables into two independent normally distributed random variables.

### Quantization error

In analog-to-digital conversion a quantization error occurs. This error is either due to rounding or truncation. When the original signal is much larger than one least significant bit (LSB), the quantization error is not significantly correlated with the signal, and has an approximately uniform distribution. The RMS error therefore follows from the variance of this distribution.

## Computational methods

### Sampling from a uniform distribution

There are many applications in which it is useful to run simulation experiments. Many programming languages come with implementations to generate pseudo-random numbers which are effectively distributed according to the standard uniform distribution.

If u is a value sampled from the standard uniform distribution, then the value a + (ba)u follows the uniform distribution parametrised by a and b, as described above.

## History

While the historical origins in the conception of uniform distribution are inconclusive, it is speculated that the term 'uniform' arose from the concept of equiprobability in dice games (note that the dice games would have discrete and not continuous uniform sample space). Equiprobability was mentioned in Gerolamo Cardano's Liber de Ludo Aleae, a manual written in 16th century and detailed on advanced probability calculus in relation to dice.[10]

## References

1. Dekking, Michel (2005). A modern introduction to probability and statistics : understanding why and how. London, UK: Springer. pp. 60–61. ISBN 978-1-85233-896-1.
2. Walpole, Ronald; et al. (2012). Probability & Statistics for Engineers and Scientists. Boston, USA: Prentice Hall. pp. 171–172. ISBN 978-0-321-62911-1. Explicit use of et al. in: |last= (help)
3. Park, Sung Y.; Bera, Anil K. (2009). "Maximum entropy autoregressive conditional heteroskedasticity model". Journal of Econometrics. 150 (2): 219–230. CiteSeerX 10.1.1.511.9750. doi:10.1016/j.jeconom.2008.12.014.
4. "Uniform Distribution (Continuous)". MathWorks. 2019. Retrieved November 22, 2019.
5. Illowsky, Barbara; et al. (2013). Introductory Statistics. Rice University, Houston, Texas, USA: OpenStax College. pp. 296–304. ISBN 978-1-938168-20-8. Explicit use of et al. in: |last= (help)
6. Casella & Berger 2001, p. 626
7. https://galton.uchicago.edu/~wichura/Stat304/Handouts/L18.cumulants.pdf
8. Nechval KN, Nechval NA, Vasermanis EK, Makeev VY (2002) Constructing shortest-length confidence intervals. Transport and Telecommunication 3 (1) 95-103
9. Wanke, Peter (2008). "The uniform distribution as a first practical approach to new product inventory management". International Journal of Production Economics. 114 (2): 811–819. doi:10.1016/j.ijpe.2008.04.004 via Research Gate.
10. Bellhouse, David (May 2005). "Decoding Cardano's Liber de Ludo". Historia Mathematica. 32: 180–202. doi:10.1016/j.hm.2004.04.001.
• Casella, George; Roger L. Berger (2001), Statistical Inference (2nd ed.), ISBN 978-0-534-24312-8, LCCN 2001025794


text:
# library.roots Submodule¶

## Module Summary¶

Interfaces for the NAG Mark 27.3 roots Chapter.

roots - Roots of One or More Transcendental Equations

This module is concerned with the calculation of zeros of continuous functions of one or more variables. The majority of problems considered are for real-valued functions of real variables, in which case complex equations must be expressed in terms of the equivalent larger system of real equations.

## Background¶

The module divides naturally into two parts.

### A Single Equation¶

The first deals with the real zeros of a real function of a single variable .

There are three functions with simple calling sequences. The first assumes that you can determine an initial interval within which the desired zero lies, (that is, where ), and outside which all other zeros lie. The function then systematically subdivides the interval to produce a final interval containing the zero. This final interval has a length bounded by your specified error requirements; the end of the interval where the function has smallest magnitude is returned as the zero. This function is guaranteed to converge to a simple zero of the function. (Here we define a simple zero as a zero corresponding to a sign-change of the function; none of the available functions are capable of making any finer distinction.) However, as with the other functions described below, a non-simple zero might be determined and it is left to you to check for this. The algorithm used is due to Brent (1973).

The two other functions are both designed for the case where you are unable to specify an interval containing the simple zero. One starts from an initial point and performs a search for an interval containing a simple zero. If such an interval is computed then the method described above is used next to determine the zero accurately. The other method uses a ‘continuation’ method based on a secant iteration. A sequence of subproblems is solved; the first of these is trivial and the last is the actual problem of finding a zero of . The intermediate problems employ the solutions of earlier problems to provide initial guesses for the secant iterations used to calculate their solutions.

Three other functions are also supplied. They employ reverse communication and use the same core algorithms as the functions described above.

Finally, two functions are provided to return values of Lambert’s function (sometimes known as the ‘product log’ or ‘Omega’ function), which is the inverse function of

that is, if Lambert’s function for , then is a zero of the function . One function uses the iterative method described in Barry et al. (1995) to return values from the real branches of (restricting ). The second function enforces no such restriction, and uses the approach described in Corless et al. (1996).

### Systems of Equations¶

The functions in the second part of this module are designed to solve a set of nonlinear equations in unknowns

where stands for transpose.

It is assumed that the functions are continuous and differentiable so that the matrix of first partial derivatives of the functions, the Jacobian matrix evaluated at the point , exists, though it may not be possible to calculate it directly.

The functions must be independent, otherwise there will be an infinity of solutions and the methods will fail. However, even when the functions are independent the solutions may not be unique. Since the methods are iterative, an initial guess at the solution has to be supplied, and the solution located will usually be the one closest to this initial guess.

## Available Routines¶

### Zeros of Functions of One Variable¶

The functions can be divided into two classes. There are three functions (contfn_interval_rcomm(), contfn_cntin_rcomm() and contfn_brent_rcomm()) all written in reverse communication form and three (contfn_brent_interval(), contfn_cntin() and contfn_brent()) written in direct communication form (see Direct and Reverse Communication Routines for a description of the difference between these two conventions). The direct communication functions are designed for inexperienced users and, in particular, for solving problems where the function whose zero is to be calculated, can be coded as a callback function. These functions find the zero by using the same core algorithms as the reverse communication functions. Experienced users are recommended to use the reverse communication functions directly as they permit you more control of the calculation. Indeed, if the zero-finding process is embedded in a much larger program then the reverse communication functions should always be used.

The recommendation as to which function should be used depends mainly on whether you can supply an interval containing the zero; that is, where . If the interval can be supplied, then contfn_brent() (or, in reverse communication, contfn_brent_rcomm()) should be used, in general. This recommendation should be qualified in the case when the only interval which can be supplied is very long relative to your error requirements and you can also supply a good approximation to the zero. In this case contfn_cntin() (or, in reverse communication, contfn_cntin_rcomm()) may prove more efficient (though these latter functions will not provide the error bound available from contfn_brent_rcomm()).

If an interval containing the zero cannot be supplied then you must choose between contfn_brent_interval() (or, in reverse communication, contfn_interval_rcomm() followed by contfn_brent_rcomm()) and contfn_cntin() (or, in reverse communication, contfn_cntin_rcomm()). contfn_brent_interval() first determines an interval containing the zero, and then proceeds as in contfn_brent(); it is particularly recommended when you do not have a good initial approximation to the zero. If a good initial approximation to the zero is available then contfn_cntin() is to be preferred. Since neither of these latter functions has guaranteed convergence to the zero, you are recommended to experiment with both in case of difficulty.

### Solution of Sets of Nonlinear Equations¶

The solution of a set of nonlinear equations

can be regarded as a special case of the problem of finding a minimum of a sum of squares

So the functions in submodule opt are relevant as well as the special nonlinear equations functions.

The functions for solving a set of nonlinear equations can also be divided into classes. There are six functions (sys_func_aa(), sys_func_easy(), sys_func_expert(), sparsys_func_easy(), sys_deriv_easy() and sys_deriv_expert()) all written in direct communication form and three (sys_func_aa_rcomm(), sys_func_rcomm() and sys_deriv_rcomm()) written in reverse communication form. The direct communication functions are designed for inexperienced users and, in particular, these functions require the (and possibly their derivatives) to be calculated in functions. These should be set up carefully so the NAG Library functions can work as efficiently as possible. Experienced users are recommended to use the reverse communication functions as they permit you more control of the calculation. Indeed, if the zero-finding process is embedded in a much larger program then the reverse communication functions should always be used.

The main decision you have to make is whether to supply the derivatives . It is advisable to do so if possible, since the results obtained by algorithms which use derivatives are generally more reliable than those obtained by algorithms which do not use derivatives.

sys_deriv_easy(), sys_deriv_expert() and sys_deriv_rcomm() require you to provide the derivatives, whilst sys_func_aa(), sys_func_aa_rcomm(), sys_func_easy(), sys_func_expert(), sys_func_rcomm() and sparsys_func_easy() do not. sys_func_easy(), sys_deriv_easy() and sparsys_func_easy() are easy-to-use functions; greater flexibility may be obtained using sys_func_expert() and sys_deriv_expert() (or, in reverse communication, sys_func_rcomm() and sys_deriv_rcomm()), but these have longer argument lists. sys_func_easy(), sys_func_expert(), sys_func_rcomm() and sparsys_func_easy() approximate the derivatives internally using finite differences. sys_func_aa() and sys_func_aa_rcomm() do not use derivatives at all, and may be useful when the cost of evaluating is high.

sparsys_func_easy() is an easy-to-use function specially adapted for sparse problems, that is, problems where each function depends on a small subset of the variables so that the Jacobian matrix has many zeros. It employs sparse linear algebra methods and consequently is expected to take significantly less time to complete than the other functions, especially if is large.

sys_deriv_check() is provided for use in conjunction with sys_deriv_easy(), sys_deriv_expert() and sys_deriv_rcomm() to check the user-supplied derivatives for consistency with the functions themselves. You are strongly advised to make use of this function whenever sys_deriv_easy(), sys_deriv_expert() or sys_deriv_rcomm() is used.

Firstly, the calculation of the functions and their derivatives should be ordered so that cancellation errors are avoided. This is particularly important in a function that uses these quantities to build up estimates of higher derivatives.

Secondly, scaling of the variables has a considerable effect on the efficiency of a function. The problem should be designed so that the elements of are of similar magnitude. The same comment applies to the functions, i.e., all the should be of comparable size.

The accuracy is usually determined by the accuracy arguments of the functions, but the following points may be useful.

1. Greater accuracy in the solution may be requested by choosing smaller input values for the accuracy arguments. However, if unreasonable accuracy is demanded, rounding errors may become important and cause a failure.

2. Some idea of the accuracies of the may be obtained by monitoring the progress of the function to see how many figures remain unchanged during the last few iterations.

3. An approximation to the error in the solution is given by where is the solution to the set of linear equations

where .

Note that the decomposition of is available from sys_func_expert() and sys_deriv_expert() (or, in reverse communication, sys_func_rcomm() and sys_deriv_rcomm()) so that

and is also provided by these functions.

4. If the functions are changed by small amounts , for , then the corresponding change in the solution is given approximately by , where is the solution of the set of linear equations

Thus one can estimate the sensitivity of to any uncertainties in the specification of , for . As noted above, the sophisticated functions sys_func_expert() and sys_deriv_expert() (or, in reverse communication, sys_func_rcomm() and sys_deriv_rcomm()) provide the decomposition of .

### Values of Lambert’s W function¶

If you require purely-real values of , these will be evaluated marginally more efficiently by lambertw_real() than lambertw_complex() owing to the differing iterative procedures used by each function.

## Tutorial¶

### Zeros of Functions of One Variable¶

Consider the problem of finding such that ; i.e., the value of that gives when .

>>> from math import exp
>>> f = lambda x: x - exp(-x)

All of the NAG zero-finding routines require you to have some basic knowledge of the structure of the function you are addressing in order that the solver can know where to begin. This knowledge must either be in the form of a rough initial guess for the zero, or as an interval on which changes sign (and, therefore, passes through zero).

In the worst case a sign-change interval is not known and the initial guess for the zero may be fairly poor:

>>> x = 1.

If so, use contfn_brent_interval():

>>> from naginterfaces.library import roots
>>> soln = roots.contfn_brent_interval(x, f)

Upon successful exit contains the computed zero and the bounding interval that was used internally:

>>> print('Final value of x is {:1.8f}.'.format(soln.x))
Final value of x is 0.56714329.
>>> print(
...     'Search interval was [' +
...     ', '.join(['{:1.8f}']*2).format(soln.a, soln.b) +
...     '].'
... )
Search interval was [0.50000000, 0.90000000].

For an alternative approach that uses continuation, more suitable for when a better initial guess is known, use contfn_cntin():

>>> x = 0.6
>>> soln = roots.contfn_cntin(x, f)

The continuation algorithm returns only the final :

>>> print('Final value of x is {:1.8f}.'.format(soln))
Final value of x is 0.56714329.

### Solution of Sets of Nonlinear Equations¶

Consider the following tridiagonal system of nine equations in nine variables:

as a vector function

>>> def f(x):
...     fvec = (3. - 2.*x)*x + 1.
...     fvec[1:] = fvec[1:] - x[:n-1]
...     fvec[:n-1] = fvec[:n-1] - 2.*x[1:]
...     return fvec

for which we seek such that .

As with scalar root-finding above, the NAG solver needs you to supply a rough initial guess for .

>>> import numpy as np
>>> n = 9; x = -1.*np.ones(n)

To solve a system of nonlinear equations, again the NAG routines can benefit from your knowledge of the structure of the problem as outlined above.

If you are happy for the Jacobian of the system to be approximated internally, use, say, sys_func_easy():

>>> from naginterfaces.library import roots
>>> soln = roots.sys_func_easy(f, x)
>>> def show_soln():
...     print(
...         'Final approximate solution is\n'
...         '(\n' +
...         '\n'.join(['  {:.4f},']*n).format(*soln.x) +
...         '\n)'
...     )
>>> show_soln()
Final approximate solution is
(
-0.5707,
-0.6816,
-0.7017,
-0.7042,
-0.7014,
-0.6919,
-0.6658,
-0.5960,
-0.4164,
)

The exact Jacobian of the system in question can be represented using the following Python function:

>>> import numpy as np
>>> def fjac(x):
...     fj = np.zeros((len(x), len(x)))
...     fj[0, 0] = 3. - 4.*x[0]
...     fj[0, 1] = -2.
...     for i in range(1, n-2):
...         fj[i, i-1] = -1.
...         fj[i, i] = 3. - 4.*x[i]
...         fj[i, i+1] = -2.
...     fj[n-1, n-2] = -1.
...     fj[n-1, n-1] = 3. - 4.*x[n-1]
...     return fj

Function sys_deriv_easy() can be used to take advantage of these exact derivatives:

>>> from naginterfaces.library import roots
>>> soln = roots.sys_deriv_easy(f, fjac, x)
>>> show_soln()
Final approximate solution is
(
-0.5707,
-0.6816,
-0.7017,
-0.7042,
-0.7014,
-0.6919,
-0.6658,
-0.5960,
-0.4164,
)

When the Jacobian is known to contain a large number of zero entries, but isn’t available in an exact form, sparsys_func_easy() may be appropriate. With this solver we are asked to compute values of the system for a supplied base- index mask. (In the contrived example below we use precomputed ‘full’ values for the system, using our definition of above, to populate the requested entries of the return vector. For a real-world sparse problem where memory usage is a concern this would clearly be counterproductive.) The function vector is updated in place in the callback, which further optimizes memory usage.

>>> def sparsf(indf, x, fvec):
...     fv = f(x)
...     for ind in indf:
...        fvec[ind-1] = fv[ind-1]

For our system the number of nonzeros is . Communicating this in the NAG call can reduce memory requirements.

>>> from naginterfaces.library import roots
>>> soln = roots.sparsys_func_easy(sparsf, x, nnz=3*n-2)
>>> show_soln()
Final approximate solution is
(
-0.5707,
-0.6816,
-0.7017,
-0.7042,
-0.7014,
-0.6919,
-0.6658,
-0.5960,
-0.4164,
)

Computational efficiency may also be enhanced by using the fixed-point iteration routine sys_func_aa() and its Anderson acceleration facility. The potential downsides are that convergence is less reliable than with the Powell hybrid method used in the solvers shown previously, and that the system may need to be rearranged into the form for best performance:

>>> import numpy as np
>>> def fixedf(x):
...     fvec = np.zeros(n)
...     fvec[0] = (2.*x[1] - 1.)/(3. - 2.*x[0])
...     for i in range(1, n-1):
...         fvec[i] = (2.*x[i+1] - 1. + x[i-1])/(3. - 2.*x[i])
...     fvec[n-1] = (x[n-2] - 1.)/(3. - 2.*x[n-1])
...     return fvec - x

The size of the iteration subspace cache to use is supplied via argument :

>>> from naginterfaces.library import roots
>>> soln = roots.sys_func_aa(fixedf, x, m=4)
>>> show_soln()
Final approximate solution is
(
-0.5707,
-0.6816,
-0.7017,
-0.7042,
-0.7014,
-0.6919,
-0.6658,
-0.5960,
-0.4164,
)

naginterfaces.library.examples.roots :

This subpackage contains examples for the roots module. See also the Examples subsection.

## Functionality Index¶

Lambert’s function

Zeros of functions of one variable

direct communication

binary search followed by Brent algorithm: contfn_brent_interval()

Brent algorithm: contfn_brent()

continuation method: contfn_cntin()

reverse communication

Zeros of functions of several variables

checking function

checks user-supplied Jacobian: sys_deriv_check()

direct communication

Anderson acceleration: sys_func_aa()

easy-to-use

derivatives required: sys_deriv_easy()

no derivatives required: sys_func_easy()

no derivatives required, sparse: sparsys_func_easy()

sophisticated

derivatives required: sys_deriv_expert()

no derivatives required: sys_func_expert()

reverse communication

Anderson acceleration: sys_func_aa_rcomm()

sophisticated

derivatives required: sys_deriv_rcomm()

no derivatives required: sys_func_rcomm()

For full information please refer to the NAG Library document

https://www.nag.com/numeric/nl/nagdoc_27.3/flhtml/c05/c05intro.html

## Examples¶

naginterfaces.library.examples.roots.contfn_brent_interval_ex.main()[source]

Finds a zero of a continuous function using Brent’s algorithm.

>>> main()
naginterfaces.library.roots.contfn_brent_interval Python Example Results.
Implied volatility for an option with a given target price.
The implied volatility for the target price of 12.35008695 is 0.09000000.

naginterfaces.library.examples.roots.sys_func_rcomm_ex.main()[source]

Find a solution of a system of nonlinear equations.

Demonstrates promoting a NagAlgorithmicWarning to an exception.

>>> main()
naginterfaces.library.roots.sys_func_rcomm Python Example Results.
Solution of a system of nonlinear equations.
Final 2-norm of the residuals after 11 iterations is 1.19264e-08.
Final approximate solution is
(
-0.5707,
-0.6816,
-0.7017,
-0.7042,
-0.7014,
-0.6919,
-0.6658,
-0.5960,
-0.4164,
)