1+ using System ;
2+ using System . Collections . Generic ;
3+
4+ using Meta . Numerics ;
5+ using Meta . Numerics . Analysis ;
6+ using Meta . Numerics . Matrices ;
7+
8+ namespace Examples {
9+
10+ public static class Analysis {
11+
12+ [ ExampleMethod ]
13+ public static void Integration ( ) {
14+
15+ // This is the area of the upper half-circle, so the value should be pi/2.
16+ IntegrationResult i = FunctionMath . Integrate ( x => Math . Sqrt ( 1.0 - x * x ) , - 1.0 , + 1.0 ) ;
17+ Console . WriteLine ( $ "i = { i . Estimate } after { i . EvaluationCount } evaluations.") ;
18+
19+ Interval zeroToPi = Interval . FromEndpoints ( 0.0 , Math . PI ) ;
20+ IntegrationResult watson = MultiFunctionMath . Integrate (
21+ x => 1.0 / ( 1.0 - Math . Cos ( x [ 0 ] ) * Math . Cos ( x [ 1 ] ) * Math . Cos ( x [ 2 ] ) ) ,
22+ new Interval [ ] { zeroToPi , zeroToPi , zeroToPi }
23+ ) ;
24+
25+ IntegrationResult i2 = FunctionMath . Integrate (
26+ x => Math . Sqrt ( 1.0 - x * x ) , - 1.0 , + 1.0 ,
27+ new IntegrationSettings ( ) { RelativePrecision = 1.0E-4 }
28+ ) ;
29+
30+ // This integrand has a log singularity at x = 0.
31+ // The value of this integral is the Catalan constant.
32+ IntegrationResult soft = FunctionMath . Integrate (
33+ x => - Math . Log ( x ) / ( 1.0 + x * x ) , 0.0 , 1.0
34+ ) ;
35+
36+ // This integral has a power law singularity at x = 0.
37+ // The value of this integral is 4.
38+ IntegrationResult hard = FunctionMath . Integrate (
39+ x => Math . Pow ( x , - 3.0 / 4.0 ) , 0.0 , 1.0 ,
40+ new IntegrationSettings ( ) { RelativePrecision = 1.0E-6 }
41+ ) ;
42+
43+ // The value of this infinite integral is sqrt(pi)
44+ IntegrationResult infinite = FunctionMath . Integrate (
45+ x => Math . Exp ( - x * x ) , Double . NegativeInfinity , Double . PositiveInfinity
46+ ) ;
47+
48+ Func < IReadOnlyList < double > , double > distance = z => {
49+ ColumnVector x = new ColumnVector ( z [ 0 ] , z [ 1 ] , z [ 2 ] ) ;
50+ ColumnVector y = new ColumnVector ( z [ 3 ] , z [ 4 ] , z [ 5 ] ) ;
51+ ColumnVector d = x - y ;
52+ return ( d . Norm ( ) ) ;
53+ } ;
54+
55+ Interval oneBox = Interval . FromEndpoints ( 0.0 , 1.0 ) ;
56+ Interval [ ] sixBox = new Interval [ ] { oneBox , oneBox , oneBox , oneBox , oneBox , oneBox } ;
57+
58+ IntegrationSettings settings = new IntegrationSettings ( ) {
59+ RelativePrecision = 1.0E-4 ,
60+ AbsolutePrecision = 0.0 ,
61+ Listener = r => {
62+ Console . WriteLine ( $ "Estimate { r . Estimate } after { r . EvaluationCount } evaluations.") ;
63+ }
64+ } ;
65+
66+ IntegrationResult numeric = MultiFunctionMath . Integrate ( distance , sixBox , settings ) ;
67+ Console . WriteLine ( $ "The numeric result is { numeric . Estimate } .") ;
68+
69+ double analytic = 4.0 / 105.0 + 17.0 / 105.0 * Math . Sqrt ( 2.0 ) - 2.0 / 35.0 * Math . Sqrt ( 3.0 )
70+ + Math . Log ( 1.0 + Math . Sqrt ( 2.0 ) ) / 5.0 + 2.0 / 5.0 * Math . Log ( 2.0 + Math . Sqrt ( 3.0 ) )
71+ - Math . PI / 15.0 ;
72+ Console . WriteLine ( $ "The analytic result is { analytic } .") ;
73+
74+ }
75+
76+ [ ExampleMethod ]
77+ public static void Optimization ( ) {
78+
79+ Interval range = Interval . FromEndpoints ( 0.0 , 6.28 ) ;
80+ Extremum min = FunctionMath . FindMinimum ( x => Math . Sin ( x ) , range ) ;
81+ Console . WriteLine ( $ "Found minimum of { min . Value } at { min . Location } .") ;
82+ Console . WriteLine ( $ "Required { min . EvaluationCount } evaluations.") ;
83+
84+ MultiExtremum rosenbrock = MultiFunctionMath . FindLocalMinimum (
85+ x => MoreMath . Sqr ( 2.0 - x [ 0 ] ) + 100.0 * MoreMath . Sqr ( x [ 1 ] - x [ 0 ] * x [ 0 ] ) ,
86+ new ColumnVector ( 0.0 , 0.0 )
87+ ) ;
88+ ColumnVector xm = rosenbrock . Location ;
89+ Console . WriteLine ( $ "Found minimum of { rosenbrock . Value } at ({ xm [ 0 ] } ,{ xm [ 1 ] } ).") ;
90+ Console . WriteLine ( $ "Required { rosenbrock . EvaluationCount } evaluations.") ;
91+
92+ Func < IReadOnlyList < double > , double > leviFunction = z => {
93+ double x = z [ 0 ] ;
94+ double y = z [ 1 ] ;
95+ return (
96+ MoreMath . Sqr ( MoreMath . Sin ( 3.0 * Math . PI * x ) ) +
97+ MoreMath . Sqr ( x - 1.0 ) * ( 1.0 + MoreMath . Sqr ( MoreMath . Sin ( 3.0 * Math . PI * y ) ) ) +
98+ MoreMath . Sqr ( y - 1.0 ) * ( 1.0 + MoreMath . Sqr ( MoreMath . Sin ( 2.0 * Math . PI * y ) ) )
99+ ) ;
100+ } ;
101+
102+ Interval [ ] leviRegion = new Interval [ ] {
103+ Interval . FromEndpoints ( - 10.0 , + 10.0 ) ,
104+ Interval . FromEndpoints ( - 10.0 , + 10.0 )
105+ } ;
106+
107+ MultiExtremum levi = MultiFunctionMath . FindGlobalMinimum ( leviFunction , leviRegion ) ;
108+
109+ ColumnVector zm = levi . Location ;
110+ Console . WriteLine ( $ "Found minimum of { levi . Value } at ({ zm [ 0 ] } ,{ zm [ 1 ] } ).") ;
111+ Console . WriteLine ( $ "Required { levi . EvaluationCount } evaluations.") ;
112+
113+ // Select a dimension
114+ int n = 6 ;
115+
116+ // Define a function that takes 2n polar coordinates in the form
117+ // phi_1, theta_1, phi_2, theta_2, ..., phi_n, theta_n, computes
118+ // the sum of the potential energy 1/d for all pairs.
119+ Func < IReadOnlyList < double > , double > thompson = ( IReadOnlyList < double > v ) => {
120+ double e = 0.0 ;
121+ // iterate over all distinct pairs of points
122+ for ( int i = 0 ; i < n ; i ++ ) {
123+ for ( int j = 0 ; j < i ; j ++ ) {
124+ // compute the chord length between points i and j
125+ double dx = Math . Cos ( v [ 2 * j ] ) * Math . Cos ( v [ 2 * j + 1 ] ) - Math . Cos ( v [ 2 * i ] ) * Math . Cos ( v [ 2 * i + 1 ] ) ;
126+ double dy = Math . Cos ( v [ 2 * j ] ) * Math . Sin ( v [ 2 * j + 1 ] ) - Math . Cos ( v [ 2 * i ] ) * Math . Sin ( v [ 2 * i + 1 ] ) ;
127+ double dz = Math . Sin ( v [ 2 * j ] ) - Math . Sin ( v [ 2 * i ] ) ;
128+ double d = Math . Sqrt ( dx * dx + dy * dy + dz * dz ) ;
129+ e += 1.0 / d ;
130+ }
131+ }
132+ return ( e ) ;
133+ } ;
134+
135+ // Limit all polar coordinates to their standard ranges.
136+ Interval [ ] box = new Interval [ 2 * n ] ;
137+ for ( int i = 0 ; i < n ; i ++ ) {
138+ box [ 2 * i ] = Interval . FromEndpoints ( - Math . PI , Math . PI ) ;
139+ box [ 2 * i + 1 ] = Interval . FromEndpoints ( - Math . PI / 2.0 , Math . PI / 2.0 ) ;
140+ }
141+
142+ // Use settings to monitor proress toward a rough solution.
143+ MultiExtremumSettings roughSettings = new MultiExtremumSettings ( ) {
144+ RelativePrecision = 0.05 , AbsolutePrecision = 0.0 ,
145+ Listener = r => {
146+ Console . WriteLine ( $ "Minimum { r . Value } after { r . EvaluationCount } evaluations.") ;
147+ }
148+ } ;
149+ MultiExtremum roughThompson = MultiFunctionMath . FindGlobalMinimum ( thompson , box , roughSettings ) ;
150+
151+ // Use settings to monitor proress toward a refined solution.
152+ MultiExtremumSettings refinedSettings = new MultiExtremumSettings ( ) {
153+ RelativePrecision = 1.0E-5 , AbsolutePrecision = 0.0 ,
154+ Listener = r => {
155+ Console . WriteLine ( $ "Minimum { r . Value } after { r . EvaluationCount } evaluations.") ;
156+ }
157+ } ;
158+ MultiExtremum refinedThompson = MultiFunctionMath . FindLocalMinimum ( thompson , roughThompson . Location , refinedSettings ) ;
159+
160+ Console . WriteLine ( $ "Minimum potential energy { refinedThompson . Value } .") ;
161+ Console . WriteLine ( $ "Required { roughThompson . EvaluationCount } + { refinedThompson . EvaluationCount } evaluations.") ;
162+
163+ /*
164+ // Define a function that takes 2n coordinates x1, y1, x2, y2, ... xn, yn
165+ // and finds the smallest distance between two coordinate pairs.
166+ Func<IReadOnlyList<double>, double> function = (IReadOnlyList<double> x) => {
167+ double sMin = Double.MaxValue;
168+ for (int i = 0; i < n; i++) {
169+ for (int j = 0; j < i; j++) {
170+ double s = MoreMath.Hypot(x[2 * i] - x[2 * j], x[2 * i + 1] - x[2 * j + 1]);
171+ if (s < sMin) sMin = s;
172+ }
173+ }
174+ return (sMin);
175+ };
176+
177+ // Limit all coordinates to the unit box.
178+ Interval[] box = new Interval[2 * n];
179+ for (int i = 0; i < box.Length; i++) box[i] = Interval.FromEndpoints(0.0, 1.0);
180+
181+ // Use settings to monitor proress toward a rough solution.
182+ MultiExtremumSettings roughSettings = new MultiExtremumSettings() {
183+ RelativePrecision = 1.0E-2, AbsolutePrecision = 0.0,
184+ Listener = r => {
185+ Console.WriteLine($"Minimum {r.Value} after {r.EvaluationCount} evaluations.");
186+ }
187+ };
188+ MultiExtremum roughMaximum = MultiFunctionMath.FindGlobalMaximum(function, box, roughSettings);
189+
190+ // Use settings to monitor proress toward a rough solution.
191+ MultiExtremumSettings refinedSettings = new MultiExtremumSettings() {
192+ RelativePrecision = 1.0E-8, AbsolutePrecision = 0.0,
193+ Listener = r => {
194+ Console.WriteLine($"Minimum {r.Value} after {r.EvaluationCount} evaluations.");
195+ }
196+ };
197+ MultiExtremum refinedMaximum = MultiFunctionMath.FindLocalMaximum(function, roughMaximum.Location, refinedSettings);
198+ */
199+ }
200+
201+ [ ExampleMethod ]
202+ public static void IntegrateOde ( ) {
203+
204+ Func < double , double , double > rhs = ( x , y ) => - x * y ;
205+ OdeResult sln = FunctionMath . IntegrateOde ( rhs , 0.0 , 1.0 , 2.0 ) ;
206+ Console . WriteLine ( $ "Numeric solution y({ sln . X } ) = { sln . Y } .") ;
207+ Console . WriteLine ( $ "Required { sln . EvaluationCount } evaluations.") ;
208+ Console . WriteLine ( $ "Analytic solution y({ sln . X } ) = { Math . Exp ( - MoreMath . Sqr ( sln . X ) / 2.0 ) } ") ;
209+
210+ // Lotka-Volterra equations
211+ double A = 0.1 ;
212+ double B = 0.02 ;
213+ double C = 0.4 ;
214+ double D = 0.02 ;
215+ Func < double , IReadOnlyList < double > , IReadOnlyList < double > > lkRhs = ( t , y ) => {
216+ return ( new double [ ] {
217+ A * y [ 0 ] - B * y [ 0 ] * y [ 1 ] , D * y [ 0 ] * y [ 1 ] - C * y [ 1 ]
218+ } ) ;
219+ } ;
220+ MultiOdeSettings lkSettings = new MultiOdeSettings ( ) {
221+ Listener = r => { Console . WriteLine ( $ "t={ r . X } rabbits={ r . Y [ 0 ] } , foxes={ r . Y [ 1 ] } ") ; }
222+ } ;
223+ MultiFunctionMath . IntegrateOde ( lkRhs , 0.0 , new double [ ] { 20.0 , 10.0 } , 50.0 , lkSettings ) ;
224+
225+ Func < double , IReadOnlyList < double > , IReadOnlyList < double > > rhs1 = ( x , u ) => {
226+ return new double [ ] { u [ 1 ] , - u [ 0 ] } ;
227+ } ;
228+ MultiOdeSettings settings1 = new MultiOdeSettings ( ) { EvaluationBudget = 100000 } ;
229+ MultiOdeResult result1 = MultiFunctionMath . IntegrateOde (
230+ rhs1 , 0.0 , new double [ ] { 0.0 , 1.0 } , 500.0 , settings1
231+ ) ;
232+ double s1 = MoreMath . Sqr ( result1 . Y [ 0 ] ) + MoreMath . Sqr ( result1 . Y [ 1 ] ) ;
233+ Console . WriteLine ( $ "y({ result1 . X } ) = { result1 . Y [ 0 ] } , (y)^2 + (y')^2 = { s1 } ") ;
234+ Console . WriteLine ( $ "Required { result1 . EvaluationCount } evaluations.") ;
235+
236+ Func < double , double , double > rhs2 = ( x , y ) => - y ;
237+ OdeSettings settings2 = new OdeSettings ( ) { EvaluationBudget = 100000 } ;
238+ OdeResult result2 = FunctionMath . IntegrateConservativeOde (
239+ rhs2 , 0.0 , 0.0 , 1.0 , 500.0 , settings2
240+ ) ;
241+ double s2 = MoreMath . Sqr ( result2 . Y ) + MoreMath . Sqr ( result2 . YPrime ) ;
242+ Console . WriteLine ( $ "y({ result2 . X } ) = { result2 . Y } , (y)^2 + (y')^2 = { s2 } ") ;
243+ Console . WriteLine ( $ "Required { result2 . EvaluationCount } evaluations") ;
244+
245+ Console . WriteLine ( MoreMath . Sin ( 500.0 ) ) ;
246+ }
247+
248+ }
249+
250+ }
0 commit comments