Approach of the value of a rent when non-central moments of the capitalization factor are known: an R application with interest rates following normal and beta distributions Ratio Mathematica Volume 39, 2020, pp. 33-54 33 Modelling the shape of sunspot cycle using a modified Maxwell-Boltzmann probability distribution function Amaranathan Sabarinath * Girija Puthumana Beena † Ajimandiram Krishnankuttynair Anilkumar‡ Abstract The 11-year sunspot number cycle has been a fascinating phenomenon for many scientists in the last three centuries. Various mathematical functions have been used for modelling the 11-year sunspot number cycles. In this paper, we present a new model, which is derived from the well known Maxwell-Boltzmann probability distribution function. A modification has been carried out by introducing a new parameter, called area parameter to model sunspot number cycle using Maxwell- Boltzmann probability distribution function. This parameter removes the normality condition possessed by probability density function, and fits an arbitrary sunspot cycle of any magnitude. The new model has been fitted in the actual monthly averaged sunspot cycles and it is found that, the Hathaway, Wilson and Reichmann measure, the goodness of fit is high. The estimated parameters of the sunspot number cycles 1 to 24 have been presented in this paper. A Monte Carlo based simple random search is used for nonlinear parameter estimation. The Prediction has been carried out for the next sunspot number cycle 25 through a model by averaging of recent cycle's model parameters. This prediction can be used for simulating a more realistic sunspot cycle profile. Through extensive Monte Carlo simulations, a large number of sunspot cycle profiles could be generated and these can be used in the studies of the orbital dynamics. Keywords: solar cycle; modelling; sunspot number 2010 AMS subject classification: 70F15; 97M10.§ * Vikram Sarabhai Space Centre, Thiruvananthapuram, Kerala, India; a_sabarinath@yahoo.co.in. † Assistant Professor, St.Gregorios College, Kottarakkara, Kerala, India; beenamabhi@gmail.com. ‡ Director, DSSAM, ISRO Head Quarters, Bangalore, Karnataka, India; ak_anilkumar@isro.gov.in. Sabarinath.A, Beena.G.P, Anilkumar.A.K 34 1 Introduction To know in advance the multitude of atmospheric processes that cause concern to mankind, in particular phenomena occurring in the solar plasma receive great consideration from the scientific world. Since the 18th century, scientists are conducting systematic research on a multitude of processes caused by solar activity. Solar Activity forecasting is crucial in scientific and technological fields such as spacecraft orbital life time prediction, airline communications and geophysical applications, mainly it is the energy source behind all phenomena driving space weather. The low Earth orbiting satellites are also influenced by solar activity (Seeds,M.A,Backman,D,[2015]; Hathaway,D.H., [2010]). However, predicting the solar cycle is challenging on the basis of time series of various proposed indicators, due to the high frequency contents, noise contamination, high dispersion level and high variability both in phase and amplitude. The prediction of solar activity is complicated by the lack of a quantitative theoretical model of the sun's magnetic cycle. The effect of solar activity is greater on space activities especially on the operations of low Earth orbiting satellites which provide significant contribution in communication, national defence and Earth mapping. Such satellites also handle a large quantity of scientific data. During higher solar activity, the maximum ultraviolet rays are emitted from the sun that heat up Earth's upper atmosphere, and expands the atmosphere. This affects the life time of operational space crafts in the low earth orbits (Whitlock,D, [2006]). Therefore better predictions of solar activity are essential to help spacecraft mission planning and design. 2 Satellite life time estimation and re-entry prediction In spacecraft mission design, orbital life time estimation is a critical activity (Whitlock, D, [2006]). Many uncertain parameters need to be considered while doing orbital life time estimation. The upper atmospheric density variation is the primary factor which is so difficult to predict. Many studies have been taken place to model the atmospheric density accurately. Orbital life time estimation community has always been looking up for better models of atmospheric density. Atmospheric models generally use parameters such as ap or Kp, and F10.7. Solar flux receives a lot of attention because it is an important parameter in determining atmospheric density. Most predictions rely § Received on October 24th, 2020. Accepted on December 17th, 2020. Published on December 31st, 2020. doi: 10.23755/rm.v39i0.532. ISSN: 1592-7415. eISSN: 2282-8214. ©Sabarinath et al. This paper is published under the CC-BY licence agreement. Modelling the shape of sunspot cycle using a modified Maxwell- Boltzmann probability distribution function 35 on the sunspot activity happening in the sun. This has been monitored since the 17th century regularly. An empirical relationship exists between the sunspot number, R, averaged over a month, and F10.7 (David A.Vallado et.al ,[2014]). F10.7 = 63.7 + 0.728 R + 0.000 89 R2, (1) From the above equation, we can see that 10.7 cm radio flux has a base level of about 63.7 solar flux units. To understand and estimate the radio emissions effectively we can use the following equation (David A.Vallado et.al, [2014]) F10.7 = 145 + 75 COS (0.001696 t + 0.35 SIN (0.00001695)), (2) where t is the number of days from January 1, 1981. We can summarise it as, atmospheric density is directly related to the solar flux, which in turn can be related to the solar activity. Studies done by different scientists and academicians shows that solar activity and solar flux have affirmed relation, a monthly estimate of F10.7 and sunspot number has been well established. Predicting the solar flux accurately can generate more accurate atmospheric density models that will help in fine tuning the fuel budget for longer satellite life. The discussion went so far reminds that the accurate prediction of the life time requires a very good predicted solar flux profile. In turn, it is sufficient to have a predicted sunspot number cycle. Since, via equation (1) one can transform sunspot numbers into solar flux. In this paper we try to predict sunspot number cycle in a simple and powerful technique. Initially, we model the sunspot cycle using a skew-symmetric probability distribution. The Maxwell-Boltzmann distribution is considered for this purpose. Then a preliminary level prediction is proposed as an average (mean) cycle of some recent cycles. Then a varying error band is derived from the past cycles. Within this error profiles, via Monte Carlo sampling, the predicted averaged cycle is transformed into many profiles. Sample profiles are taken and plotted. Before venturing into the details, a brief review of sunspot data and review some of the recent models are provided. 3 Sunspot number cycles and sunspot number data In 1848 the Swiss astronomer Johann Rudolph Wolff introduced actual measurements of sunspot number. His method uses still today. Total number Sabarinath.A, Beena.G.P, Anilkumar.A.K 36 of spots visible on the face of the sun is 'n' and the number of groups into which they cluster is 'g' then the sunspot number Rn is defined as Rn = 10 g + n. (3) To compensate the observational limitations like Earth's atmosphere variability above the observing site and sun's rotation, each daily sunspot number is computed as a weighted average of measurements made from a network of observatories. The 11-year cyclic variation in the sunspot numbers was first noted by Schwabe, M., [1844]. In 1848 Rudolf Wolf at Swiss Federal Observatory in Zurich, Switzerland devised his measure of sunspot numbers that continues to this day as the International sunspot number. Wolf recognised that it is far easier to identify sunspot groups than to identify each individual sunspot. This relative sunspot number,Rz with emphasis on sunspot groups is defined as, Rz = k (10 g + n), (4) Where k the correction factor for the observer, g is the number of identified sunspot groups, and n is the number of individual sunspots. These sunspot numbers are called the Zurich or International sunspot numbers have been obtained daily since 1848. Sunspot cycle time series is one of the longest time series which was studied by many experts for various reasons. First of all, this time series is non-stationary, cyclic and highly nonlinear in the time domain. In the present study, the prediction of sunspot cycles is carried out with the monthly averaged sunspot number values. The monthly averaged sunspot data were available from, http://www.sidc.be/silso/versionarchive at the royal observatory, Belgium is being used for the present study. It may be noted that, the scientific community recently recalibrated the entire historical sunspot number record and that SILSO (Sunspot Index and Long-term Solar Observations) maintains this new definitive record as well as the original version of sunspot numbers. Figure 1: Sunspot cycle evolution-Monthly averaged sunspot numbers from the year 1749 to December 2016. 1749 1791 1832 1874 1916 1957 1999 0 50 100 150 200 250 300 M o n th ly a v e ra g e d s u n s p o t n u m b e rs Time (Year) Modelling the shape of sunspot cycle using a modified Maxwell- Boltzmann probability distribution function 37 4 Existing models of sunspot number cycles Several mathematical functions were introduced to model the shape of the sunspot number cycle. Due to the exponential rise and decay, the exponential function was used by Nordemann, [1992], Nordemann, et.al.,[1992]. The bell shaped nature of the sunspot cycle was explored by Hathaway et.al.,[1994]. Few statistical probability distribution functions were also proposed for the shape modelling by various authors. De Mayer, F.,[1981], proposed a model using periodic functions. In prediction, averaged models are used as an initial estimate of the future cycle. We have an exhaustive list of details and voluminous data literature available at hand pertaining to the attempts to predict the future behaviour of solar activity (Hathaway, et.al., [1999]). It can be categorised under five heads, based on the nature of the prediction methods. They are: 1) Curve fitting, 2) Precursor, 3) Spectral, 4) Neural Networks and 5) Climatology (Sello, S.,[2001]). McNish-Lincoln curve fitting was the first attempt on the methodology of curve fitting (de Meyer,[1981], McNish, A.G., Lincoln, J.V.,[1949]). Over the years, various techniques and models have been proposed by several authors working in the field for the prediction of the nonlinear behaviour of sunspot cycles. The first breakthrough in the field of modelling the shape of the sunspot cycles by fitting an exponential function over the sunspot number cycle time series was due to Nordemann,[1992]. In this method, fitting the rise to maximum and the fall to minimum were fitted with a function of exponential function demanding six free parameters. Later a modified version of F-distribution density function with five parameters was proposed by Elling and Schwentek[1992]. Nordemann's[1992] method suggests exponential fitting and explain the solar behaviour. Hathaway, Wilson, and Reichmann[1994] substantiated the superiority of a new model along with a measure for the goodness of fit. Number of free parameters in this model is reduced to four. All these models introduce high amount of error in the prediction, due to the incompetence to fit the peak locations of the sunspot cycle. The continuous nature of the model at the high solar activity period contributes a large amount of uncertainty and hence in the applications such as the orbital re-entry predictions these models are not suitable. The next subsection surveys the literature pertaining to some models, especially on the shape of sunspot cycles. 4.1 Stewart and Panofsky model Stewart and Panofsky [1938] proposed a function for the shape of the cycle with the form Sabarinath.A, Beena.G.P, Anilkumar.A.K 38 R(t) = a(t − t0) be−c(t−t0), (5) where a, b, c, and t0 are parameters that vary from cycle to cycle. The important thing to be noticed is that, this model gives a power law for the rising phase of a cycle and an exponential for the declining phase of a cycle. The model parameters for cycle 1to 16 were computed and there by the maximum amplitude, the epoch of the peak sunspot number, etc. was predicted. 4.2 Nordemann model Nordemann used the solution of the differential equation dN dt = KN, in analogy with the nuclear decay process. Thus the declining phase of a sunspot cycle is represented by: N = N0e Kt K < 0 (6) and the solution of dN dt = A + KN, is used to represent the ascent phase of a sunspot cycle. Thus the model for the ascent phase is: N = A K (1 − eKt) K < 0 (7 ) Where N represents sunspot numbers, K decay constant and A a production parameter. The estimated values of the parameters N0, K and A for all the 22 sunspot cycles were given in Nordemann [1992]. 4.3 Elling and Schwentek model Instead of using yearly means, quarterly averages of sunspot numbers were utilised by Elling and Schwentek[1992] for optimal fitting of each cycle. They used a modified F-distribution density function that required five free parameters. This approach is much more worth than the previous models. In this model fitting concluded only for modern era of sunspot cycles (10 to 21). By considering the maxima and minima of mean sunspot number as a function of time, affinity can be observed in each cycles. While considering different sunspot cycles the ascending phase take dwindle time than the descending phase, that means Starting from a minimum, time taken for reaching the maximum is always shorter as compared to the time from maximum down to minimum . Modelling the shape of sunspot cycle using a modified Maxwell- Boltzmann probability distribution function 39 They explained very effectively, ascending and descending branches of the various cycle curves have curvatures which are rather similar to those of the F- distribution curves. For this reason, each sunspot cycles from cycle 10 to cycle 21 has been approximated by a modified F-distribution, f(t) which is defined by: f(t) = P4 Γ [ P2 + P3 2 ] Γ ( P2 2 ) Γ ( P3 2 ) P 2 P2 2 P 3 P3 2 [P1(t + P5)] P2 2 −1 [P3 + P2P1(t + P5)] (P2+P3) 2 , (8) where t is the time and Γ(x) is the gamma function. P1 is the length or duration of the sunspot cycle, that is, the time interval from one minimum to the next, P2 to the curvature of the ascending branch of f(t), P3 to the curvature of the descending branch of f(t), P4 to the amplitude of the maximum of f(t), P5 to the time shift of the f(t) curve. Through least square fit all the five parameters are estimated. 4.4 Hathaway, Wilson, and Reichmann model Hathaway et.al [1994] suggested a model with free parameters fewer than the models which we had come across. They utilised a four-parameter quasi- Planck function to fit the monthly mean sunspot numbers of a solar cycle, similar to that of Stewart and Panofsky[1938]. But the only difference we can see that a fixed power law for the initial rise of the sunspot cycle and the phase starting from maximum down to minimum can be well represented by a function that decreases as e−t 2 . By combining these, the model as a function of time can be written as: f(t) = a(t − t0) 3 e [ (t−t0) 2 b2 ] − c , (9) This model has four parameters. a represents the amplitude and is directly related to the rate of rise from minimum; b is related to the time in months from minimum to maximum; c gives the asymmetry of the cycle; and a starting time t0 . Along with the early detection of parameters to predict the solar activity they examine the relationship between the parameters. It is similar to the Plank function but contains four free parameters and has a more rapid decrease after maximum, but causes lack of accuracy. The estimation of these parameters was obtained through Levenberg-Marquardt methods (Press, W, H., [1992]). Sabarinath.A, Beena.G.P, Anilkumar.A.K 40 4.5 Volobuev’s one-parameter fit In 2009, Volobuev introduced a function of two-parameters and he refers to this as a one parameter fit. We can see that the parameters are correlated (r = 0.88) for all the 23 solar cycles. The correlation between the parameters provides the possibility of a one-parameter fit by neglecting the need to determine the best starting time. He showed that a one-parameter fit can also be derived from truncated dynamo models. Due to the unavoidable uncertainty of starting time goodness of fit value is not better as compared to the empirical fit. We can see that this model is also similar to that of Stewart and Panofsky [1938] proposed Pearson's type III curves by putting b =2 and modifying the growth multiplier and decay multiplier properly by introducing the new parameters Ts and Td. The empirical model used is written as: R = ( t − t0 Ts ) 2 e −( t−t0 Td ) 2 , (10) 4.6 Sabarinath and Anilkumar model Sabarinath and Anilkumar[2008] proposed a model consist of a mixture of Laplace distribution with six parameters (later reduced to two). This model fits the multiple sharp peaks in a solar cycle. The model for a generic cycle is: F = A1 33.2 exp ( −|t − 41.7| 16.6 ) + A2 46 exp ( −|t − 67.3| 23 ) , (11) where t is the time. 5 Skew symmetrical distributions Sunspot cycles are asymmetric with respect to their maxima (Hathaway, D.H., [2010]). Starting from minimum the time taken to reach maximum is 48 months and 84 months to fall back to minimum again. An average cycle can be constructed by stretching and contracting each cycle to the average length and normalising each to the average amplitude. In general, if we survey any model of the shape of the sunspot cycle, it is evident that, all functions are a product of a polynomial and a negative exponential function. Then the goodness of fit solely depends on how the model parameters are chosen in the model. In this context, we propose a skew Modelling the shape of sunspot cycle using a modified Maxwell- Boltzmann probability distribution function 41 symmetrical function from the class of skew symmetrical probability functions. 6 Maxwell-Boltzmann probability distribution function In statistical physics, Maxwell-Boltzmann distribution is a probability distribution named after the famous Scottish physicist James Clerk Maxwell and Ludwig Boltzmann. It is used in atomic physics for describing particle speeds in idealised gas. The Maxwell-Boltzmann distribution function is given as (Balakrishnan, N., Nevzorov, V.B., [2003]) f(v) = √( m 2πkT ) 3 4πv2e − mv2 2kT , (12) where m the particle mass and kT is the product of Boltzmann's constant and thermodynamic temperature. From Equation (12), if we put α = √ kT m , then the Maxwell-Boltzmann probability distribution function can be simplified as f(x; α) = 1 α3 √ 2 π x2e − x2 2α2 , (13) where the variable v is replaced with a generic random variable x with x ≥ 0 and it can be noted that the parameter α ≥ 0 is a real quantity. Typical shape of Maxwell-Boltzmann distribution is given in Figure-2, for a value of =30. One can clearly see from Figure-2 that the ascend phase is of 47 units and the descent phase is 85 units. There by, a skew symmetrical process or phenomenal could be modelled by the Maxwell-Boltzmann distribution. Our interest is in modelling the sunspot cycle. By observing all the cycles individually one can easily see that the rise time (starting minimum to maximum sunspot number) and fall time (maximum sunspot number to cycle end) are not equal or not symmetrical about the peak sunspot number occurring epoch during the 11 year sunspot cycle period. Sabarinath.A, Beena.G.P, Anilkumar.A.K 42 Figure-2. Maxwell-Boltzmann distribution for a value of =30.0 7 Modified Maxwell-Boltzman probability distribution function (MMPDF) Since equation (13) being a probability density function, we know that, mathematically the area under the probability density function is 1, that is, ∫ f(x)dx = 1, (14) ∞ −∞ So, if we want to fit this equation (13) into an arbitrary set S of N data points, S = {(xi, yi); xi ∈ R, yi ∈ R, i = 1,2, … , N}, where R is the set of real numbers, we need to de-normalise the property of f(x) given by equation (14). This is because; the area under the curve determined by the set of points in S need not be equal to one. That is, ∑ [(xi − xi−1) (yi + yi−1) 2 ] N i=2 = A, (15) where A need not be equal to 1. In this case we can modify equation (13) to fit into any arbitrary set as equation (16) by introducing a new parameter called area parameter A. f(x; α; A) = A α3 √ 2 π x2e − x2 2α2 , (16) Modelling the shape of sunspot cycle using a modified Maxwell- Boltzmann probability distribution function 43 Now, it may be noted that, ∫ f(x)dx = A, (17) ∞ −∞ Modified model for the sunspot cycles is f(x; α; A) = A α3 √ 2 π x2e − x2 2α2 , (18) where A is the area parameter. Modified Maxwell-Boltzmann distribution with a value of =30 and A=6000 is given in Figure-3. Figure-3. Modified Maxwell-Boltzmann distribution for a value of =30 and A=6000. 8 Estimation of model parameters The function in which parameters to be estimated is, f(x; α; A) = A α3 √ 2 π x2e − x2 2α2 . (19) The maximum likelihood estimate of the parameters α and A are considered to be the best unbiased, consistent and sufficient estimate of the parameters Sabarinath.A, Beena.G.P, Anilkumar.A.K 44 (Sorenson, H.W., [1980]). Practically, the least square estimate is considered to be the maximum likelihood estimate. The simple mathematical procedure to estimate the parameters is to minimise the sum of squared error function J, J = ∑ er 2 r , (20) Where er is the error. The minimum of J can be found by differentiating J with respect to the parameters α and A. In the present study, if we consider without loss of generality, a sunspot cycle having a length of 132 months( 11 year), and if we assume {sn: n = 1,2, … ,132} as the realised sunspot number values, then the J function can be written as, J = ∑[sn − f(xn, α, A)] 2 132 n=1 , (21) where, xn = 1,2, … ,132, represents the months for each n = 1,2, … . ,132. Then our objective is to compute and solve α and A from ∂J ∂α = 0, (22) ∂J ∂A = 0, (23) Analytically solving the equations (22) and (23) for α and A is not possible due to the nonlinear terms involved in the equations. Hence we go with numerical procedures for estimating the parameters. Monte Carlo based simple random search based procedure is considered here to estimate the parameters. This procedure is described below as an algorithm. Step-1. Start with a search region α and A. Let Sα and SA are the bounded search regions of α and A. Our objective is to find an α0 ∈ Sα and A0 ∈ SA, such that, Jα0,A0 = ∑[sn − f(xn, α0, A0)] 2 132 n=1 , (24) is minimum or Jα0,A0 ≤ Jα,A (25) for any α ∈ Sα and A ∈ SA. Step-2. Start with a random initial value of α in Sα and A in SA. Compute J and in each iteration keep the minimum value of J, α and A. After a very large number of iterations take the value of, α and A corresponds to the global minimum value of J. Modelling the shape of sunspot cycle using a modified Maxwell- Boltzmann probability distribution function 45 9 Fitting of MMPDF on sunspot cycles Using the method described in section 8, the model parameters are estimated for all the past 24 cycles. It is noticed that the fit is very much close to the actual sunspot numbers. This is evident in the goodness of fit computed for each of the 24 cycles, which is discussed in the next section in detail. Figure 4 and 5 shows the model and actual data of sunspot cycles 20 and 22. Figure-4. Fitting of sunspot cycle 20 by the model Sabarinath.A, Beena.G.P, Anilkumar.A.K 46 Figure-5. Fitting of sunspot cycle 22 by the model Figure-6. The parameters α and A for all the 24 cycles. Modelling the shape of sunspot cycle using a modified Maxwell- Boltzmann probability distribution function 47 10 Models of sunspot cycles 1 to 24 The estimated model for all the past 24 cycles is given in Table-1. In Figure-6, the variation trends of the parameters α and A for all the 24 cycles are plotted. It may be noted that, the average of the parameters are 36.25 units of α and 7095.76 of A. Table-1. Estimated parameters of cycles 1 to 24 Cycle No α A 1 48.76 5883.33 2 33.40 6251.65 3 30.56 7309.46 4 35.80 8619.41 5 43.79 3525.58 6 48.75 3067.09 7 48.16 5322.72 8 32.65 7552.73 9 44.70 8234.25 10 40.63 6410.82 11 33.80 7381.50 12 37.13 4433.49 13 32.55 4933.80 14 38.57 4356.00 15 36.00 5390.27 16 35.73 4882.42 17 38.95 7341.83 18 36.35 9228.41 19 33.79 11420.62 20 40.02 7959.33 21 35.49 9907.72 22 31.48 9075.22 23 38.99 8006.39 24 38.65 5023.60 Mean 1 to 24 38.11 6729.90 Mean 11 to 24 36.25 7095.76 It may be noted that variation in α is less and variation in A is more. So A is a more sensitive parameter than α. Variation in A is not much significant as its sensitivity is less. Sabarinath.A, Beena.G.P, Anilkumar.A.K 48 10.1 Goodness of fit Goodness of fit by Hathaway, Wilson, and Reichmann [1994] is measured by the following function χ = √ ( ∑ (Ri − fi) 2N i=1 si 2⁄ ) N , (26) where, Ri and si is the monthly averaged sunspot number and its standard deviation respectively , fi gives the functional fit value, N is the number of months in the cycle. Using this equation, computed χ value for all the 23 cycles. For Checking the Goodness of fit of the proposed model we have to consider other popular methods available in the literature. The second column of Table 2 gives the goodness of fit of the proposed Modified Maxwell- Boltzmann distribution function; the third and the fourth column gives the goodness of fit by three and two parameter fit of Hathaway, Wilson, and Reichmann [1994], respectively; the fifth column gives the goodness of fit by the five parameter function of Elling and Schwentek[1992] who considered cycles 10 to 21 for their study. Figure-7, shows the goodness of fit of 3 different models along with the Modified Maxwell-Boltzmann distribution function model. It may be observed from the goodness of fit value, that the present model proposed in this study has a very good fitness compared with other models. Especially the modern cycles (cycles 11 to 24) shows very good fitness for the Modified Maxwell-Boltzmann distribution function model. Figure-7. The goodness of fit of 3 different models and MMPDF model Modelling the shape of sunspot cycle using a modified Maxwell- Boltzmann probability distribution function 49 Table 2: Hathaway, Wilson and Reichmann χ -measure of the goodness of fit value computed for all the 22 sunspot cycles with different models. MMPDF shows good fit compared with other models. Cycle Number MMPDF model Three- parameter fit by Hathaway et.al. Two- parameter fit by Hathaway et.al. Elling- Schwentek F-distribution fit 1 0.69 0.71 0.75 2 1.38 1.42 1.50 3 1.64 1.70 1.56 4 0.93 0.89 0.95 5 2.87 2.34 2.50 6 1.72 1.90 2.14 7 1.80 0.94 1.01 8 1.16 0.96 0.99 9 0.86 0.99 0.97 10 0.72 0.74 0.76 0.70 11 0.75 0.88 0.83 1.35 12 2.06 2.08 2.12 2.17 13 0.70 0.90 0.91 0.90 14 0.97 1.11 1.09 1.12 15 0.80 0.88 0.89 1.16 16 0.76 0.89 0.97 0.89 17 0.98 0.86 0.87 1.10 18 1.21 1.05 1.04 1.27 19 0.90 0.91 0.89 1.61 20 0.79 0.87 0.95 0.66 21 0.94 0.89 0.89 1.11 22 0.82 1.05 1.06 23 0.79 11 Prediction of sunspot cycle 25 As an attempt to predict the sunspot cycle 25, we consider the average of the model parameters by considering cycles-11 to 24. This computed average is given in Table-1. Thus, the parameter values of cycle 25 are: α = 36.25, and A = 7095.76. Hence the model is, Sabarinath.A, Beena.G.P, Anilkumar.A.K 50 f(x; α; A) = A α3 √ 2 π x2e − x2 2α2 , (27) where, α = 36.25, and A = 7095.76. That is, f(x; 36.74; 6608.04) = 0.119 x2e−0.00038x 2 , (28) is the model for cycle-25. Figure-8 shows the shape of cycle 25 in an average sense. It may be observed that cycle 25 may peak up to 105 units and it is also fairly a slow cycle as cycle 24. Figure-8. Preliminary level prediction of sunspot cycle 25 . 12 Prediction error and simulated sunspot cycles Any prediction or forecast is partial, if it is not supplemented with a prediction error. Here, for our study we propose a prediction error band based on the statistical variation of all the cycles. For this, consider all the monthly averaged cycles. We propose the error band each month data as ±s, where s is the standard deviation of the sunspot numbers for that month. Figure-9 shows the mean along with the mean+s, the upper bound, and mean-s, the lower bound profile. Modelling the shape of sunspot cycle using a modified Maxwell- Boltzmann probability distribution function 51 Figure-9. Mean cycle from the actual monthly sunspot cycle along with the mean+s, the upper bound and mean-s, the lower bound profile Once we are having a prediction error and a prediction model, we can generate any number of forecast profile based on simple Monte Carlo method. Here we consider the envelop derived above as the envelope with 99.7% confidence or 3sigma confidence level, since all the realised cycles falls inside the proposed confidence interval band. Hence in the Monte Carlo simulation a typical profile will be generated using equation (29). sn ′ (i) = mn(i) + rand(i) × ( env(i) 3 ) , (29) where sn ′ (i) , is the simulated n-th sunspot cycle, i = 1,2, … , Cycle length, mn(i) is the model value, rand(i) is the random number and env(i) is the envelop value given in Figure-10. Sabarinath.A, Beena.G.P, Anilkumar.A.K 52 Figure-10. Simulated sunspot cycle 20 by the model The same methodology proposed in the study can be implemented to the F10.7 cm solar flux value and one can easiliy forecast an entire cycle and subsequently it can be applied in the life time computation of satellites. 13 Conclusions The 11-year sunspot number cycles have been a fascinating phenomenon for many in the last three centuries. Different mathematical models have been derived for modelling the shape of the 11-year sunspot number cycles. In the present study, we introduced a new model which is derived from the well known Maxwell-Boltzmann probability distribution function. The modification has been carried out by introducing a new parameter, called area parameter. The new model has been fitted in the original monthly averaged sunspot cycles data and it is found that a very high goodness of fit through the Hathaway, Wilson and Reichmann measure. The models estimated for all the sunspot cycles from 1 to 24 have been presented. Detailed discussion on the nonlinear parameter estimation carried out for fitting the function in the original data is also summarised. An attempt has been carried out for predicting the next sunspot cycles 25. The sunspot cycle 25 may peak up to 105 units and it is also fairly a slow cycle as the previous cycle 24. Modelling the shape of sunspot cycle using a modified Maxwell- Boltzmann probability distribution function 53 References [1].Balakrishnan, N., Nevzorov, V.B., A primer on statistical distributions, John Wiley & Sons, Inc., Hoboken, New Jersey. 2003. [2].De Mayer, F., Mathematical Modelling of the sunspot cycle, Solar Phys., Vol. 70, 259-272,1981. [3].David, A. Vallado, David Finkleman, A Critical Assessment of Satellite Drag and Atmospheric Density Modeling. 2014. [4].Elling, W., Schwentek, H., Fitting the sunspot cycles 10-21 by a modified F-distribution density function, Solar Phys., Vol. 137, 155-165, 1992. [5].Hathaway, D.H., The solar cycle, Living reviews in Solar Physics, 7, 1. 2010. [6].Hathaway, D.H., Wilson, R.M., Reichmann, E.J., J. Geophys. Res., Vol. 104, No. A10, 1999. [7].Hathaway,D.H., Wilson, R.M., Reichmann, E.J., The shape of the sunspot cycle, Solar Phys., April 1994, Vol. 151, No.1, pp 177-190.1994. [8].McNish, A.G., Lincoln, J.V., Prediction of sunspot numbers, Eos Trans. AGU, 30, 673–685.1949. [9].Nordemann, D.J.R., Sunspot number time series: exponential fitting and solar behavior, Solar Phys., Vol.141, 199-202.1992. [10].Nordemann, D.J.R., Trivedi, N.B., Sunspot number time series: exponential fitting and solar behavior, Solar Phys., Vol.142, 411-414.1992. [11].Petrovay, K., Solar cycle prediction. Living reviews in solar physics, 7(1):6.2010. [12].Press, W, H., The art of scientific computing. Cambridge university press.1992. [13].Sabarinath, A., Anilkumar, A.K., Modelling of sunspot numbers by a modified binary mixture of Laplace distribution functions, Solar Phys., Vol. 250, 183-197.2008. [14].Seeds, M.A, Backman, D, The solar system. Nelson Education.2015. [15].Sello, S., Solar cycle forecasting: a nonlinear dynamics approach. Astronomy & Astrophysics, 377(1):312–320. 2001. [16].Schwabe, H. von Herrn Hofrath Schwabe, Sonnen -Beobachtungen im jahre 1843. in dessau. Astronomische Nachrichten, 21:233. 1844. Sabarinath.A, Beena.G.P, Anilkumar.A.K 54 [17].Sorenson, H.W., Parameter Estimation-Principles and Problems, Marcel Dekker, New York.1980. [18].Stewart, J.Q., Panofsky, H.A.A., The Mathematical Characteristics of Sunspot Variations, Astrophys. J., Vol. 88, 385–407. 1938. [19].Volobuev, D.M., The shape of the sunspot cycle: A one parameter fit, Solar Phys., Vol. 258, 319-330.2009. [20].Whitlock, D, Modeling the effect of high solar activity on the orbital debris environment, Orbital Debris Quarterly News, 10(2):4.2006.