On the implementation of novel RKARMS(4,4) algorithm to study the structures of initial extrasolar giant protoplanets.

The main motivation and novel notion of this present communication is to implement the recently suggested fourth order with four stages embedded RKARMS(4,4) algorithm to examine its efficiency in reinvesting the structures of extrasolar protoplanets formed via disk instability which being presented in Paul et al. [1] (G.C. Paul, M.M. Rahman, D. Kumar, M.C. Barman, the radius spectrum of solid grains settling in gaseous giant protoplanets, Earth Sci. Inform. 6 (2013) 137-144) for the case of convective heat transfer using classical Runge-Kutta (RK) technique of order four. The results by the RKARMS(4,4) algorithm compared well with those obtained by the classical RK method of order four for any time length and found to be more suitable.

Introduction

The formation of planetary system is one of the prominent topics to humankind since the dawn of civilization and the advancement of research still demonstrates interest towards the same. The discovery of extrasolar planets has intensified interest in them, and a large volume of research works on the physical conditions prevailing inside such planets has been conducted both outside and inside our own solar system [2, 3]. There are two feasible end mechanisms in the literature for explaining the evolution process of such planets, namely (i) core accretion and (ii) disk instability [1]. In the first scenario, the formation of gas giant planets or simply gas giants starts with core formation followed by gas accretion [4, 5], whereas in the disk instability scenario, the segregation occurs prior to the core formation. This school of thought believes that under the appropriate conditions, an instability can occur in the protoplanetary disk. Such a kind of instability can then lead to the formation of self-gravitating clumps composed of dust and gas. These clumps, in turn, can contract to form gravitationally bound sub condensations [6, 7]. Though some questions arise with regard to whether the gravitational instability can be able to form stable protoplanets from gravitationally unstable gas disk, this mechanism is thought to be capable of forming giant planets rapidly in our own solar system and elsewhere [8]. It is noticeable that planetary evolutions are highly dependent on the initial configurations of the protoplanets. Unfortunately, so far, there does not exist a single model with its definite initial structure and various models in this regard report different initial profiles [9, 10, 11, 12]. It is of interest to note here that an analytic solution of the system of equations used in configuring structure of such a protoplanet is not possible [12] without any drastic assumptions. Therefore, one may need to rely on numerical techniques in this respect. It is mentionable at this juncture that planning of new algorithms always play an important role in research to search for the best solution for any real-time problems of initial valued arising in mathematical physics, which are being solved by Runge–Kutta (RK) method and by its successive developed techniques. These techniqes are applied to solve many problems in commucication and signal processing, and to analyze electronic and transistor circuits [13], etc. Shampine and Gordon [14] explored the normal order of a RK technique having the approximate number of leading terms of an infinite Taylor series that computes the trajectory of a moving point. A novel 4th order RK technique based on the root mean square formula is addressed in [15] for solving IVPs. Yaakub and Evans [16] presented a novel 4th order RK algorithm for solving ordinary differential equations of initial valued with error control (EC). Butcher [17, 18, 19] in his studies developed the best RK pair with an error estimate (ERREST). Bader [20, 21] introduced the RK–Butcher algorithm for finding out the truncation ERRESTs, intrinsic accuracies and the early recognition of stiffness in coupled differential equations arising in theoretical problems of chemistry. With the intension of overcoming step size constraints attributed by numerical stability, numerous new algorithms in this regard have been advanced in the recent past [17]. To ensure this, in [22] a novel RKARMS(4,4) algorithm, a combination of classical RK arithmetic mean method of order four (RKAM(4,4)) and the 4th order RK root mean square (RKRMS(4,4)) method, with EC in detail is introduced to solve effictively the problems in image processing under cellular neural network (CNN) model. A detailed interpretation with regard to local and global truncation errors, ERREST and control for the 4th order and four-stage RK algorithm is finally addressed in [23]. It is to be mentioned here that that an embedded technique in fact consists of two techniques built into one, the leading one is of order and the other one is that of . The difference of the field values obtained by the techniques makes available an ERREST for the leading technique. Recently, Paul et al. [24] used that RKARMS(4,4) technique to solve vertically integrated shallow water equations (SWEs) intended for numerical foreseeing of surge levels preciously along the coast of Bangladesh, where they found the method to be efficient in computing surge levels accompanying with a cyclonic storm. In exploring initial profiles of protoplanets formed via disk instability, Paul and Senthilkumar [25] used this method. They found that the method is insensitive in solving end point constraints. It is noteworthy to mention here that Paul and Senthilkumar [25] investigated the structures of initially formed protoplanets through gravitation instability assuming the conductive-radiative heat transport.

In this paper, it is our intention to test the efficiency in terms of ERREST, accuracy, solving end point contraints and computational cost for varying end points with differnet initial time steps of the novel RKARMS(4,4) algorithm in exploring the initial structures of extrasolar giant protoplanets formed via gravitational instability assuming them to be in convective equilibrium.

The reminder of the article is arranged as follows. In Section 2, the problem statement is outlined. Section 3 demonstrates the numerical procedure adopted in solving structure equations. Section 4 deals with a brief description on the embedded RKARMS(4,4) technique. Section 5 addresses the analysis of EC for the method along with local truncation error (LTE). Discussion of results with validity and effectiveness of the technique is addressed in section 6. Finally, section 7 provides conclusion, future trend and direction of our research work.

The problem statement

As in Paul et al. [1, 12], the model of the present study undertakes a spherical giant gaseous protoplanet having a solar composition with mass ranging from 0.3 to 10 M (1M = 1.8986×1030 gm). The reasoning abaft such a consideration of the range of mass is that it covers the majority of the detected mass-range of the giant extrasolar planets as well as the giant planets of our solar system [10]. Following Paul et al. [25], it is assumed that the protoplanet originated by gravitational instability is a static body, where the ideal gas law holds good and the only source of energy here is the gravitational contraction. The assumptions can be found to be supported by DeCampli and Cameron [9] and Bodenheimer et al. [26]. Now, if the energy is assumed to be transferred by the process of convection, the structure of such an object can be specified by means of the following system of equations as under [1]:

The equation expressing hydrostatic balance,

The conservation of mass equation,

The convection heat transfer equation,

The Clausius-Clapeyron equation,

In Eqs. (1), (2), (3), and (4), , and characterize the pressure, density and temperature, respectively, of the protoplanet at a distance measuring from its centre; denotes its mass interior to ; represents the ratio of the specific heats; represents the mean molecular weight; designates the mass of a hydrogen atom and denotes Boltzmann's constant; stands for the constant of universal gravitation.

Boundary conditions

Consider a sphere of infinitesimal radius at the centre of a protoplanet. Then the mass interior to can be specified bywhere can be treated as a constant, which in turn leads to set as . Eq. (5) also validates that at the surface of the protoplanets, i.e., at , . Further, initially formed protoplanets having fairly cold centre must have reasonably low surface temperature [12]. Hence, the surface temperature, in the first approximation, can be set to zero. In addition, the mass of the protoplanetary atmosphere is just a minute fraction of the total mass of the protoplanet. Thus, in the first approximation, we may take into account the surface pressure as zero [27]. Therefore, the approximated boundary conditions (BCs) can be set as

Integration of the equations

To solve a problem numerically, dimensionless groups help to scale the problem. It is to be mentioned here that non-dimensionalization is the exclusion of units fully or partially from an equation containing physical quantities by a suitable replacement of variables. This technique can make things easier and parameterize problems where measured units are involved. In this paper, as in [27], the model equations were non–dimensionalized using Schwarzschild transformations, being given by

In coordination with Eq. (4) together with the help of the above transformations, Eqs. (1), (2), and (3) can, respectively, be brought to the following forms [1, 12]:

and

While based on the specified transformations, the BCs specified by means of (6) can be set to the following form:

Now, because of the existence of singularities in Eqs. (7), (8), and (9), one cannot start the integration from either of the boundaries, the surface or centre. Thus, we need approximate boundary conditions and that can be attained through Frobenius method. As in Paul et al. [12], the following approximate surface boundary conditions are adopted:where E is a constant to be determined.

By means of these approximated values mentioned above as our initial conditions, we have solved the system of equations specified by Eqs. (7), (8), and (9) by the newly proposed RKARMS(4,4) method from a point very close to the surface downwards to a point at the close proximity of the centre. As in [12], a value to the constant E is specified and the correct value of E is estimated through calibration. In this regard, the additional BC is checked to satisfy. The best value of E is found to be E = 45.4, which is the same as that to the presented value for the parameter in [12].

In the present numerical estimate, the radii and masses of the protoplanets are supplied from the investigation due to Helled and Schubert [10]. Also, we have used , applicable for a molecular gas, while the remaining parameters involved in the problem have been considered to have their usual values.

A note on the new RK embedded RKARMS(4,4) technique

Let be the interval over which we need to find the solution to the IVPsubject to an initial condition , where the function is sufficiently differentiable in a neighborhood of the exact solution .

To solve effectively Eq. (11), a RK embedded technique can be of help [24, 25]. A general s–stage RK pair can be described by way of the extended Butcher tableau of parameters presented in Table 1.

Table 1

The Butcher tableau.

C	A
	b^T
	bˆ^T
	E^T

In the above Table 1, , and and . The values of using the two methods at x = xn+1 can be expressed, respectively, in the following forms:where represents the step size, ,. Here, the vectors , and are dimensional and the matrix is dimensional.

In the case of an embedded technique, the LTE can be calculated through . It is pertinent to note here that the LTE controls the step size h [22]. The Butcher array form for four stages method is presented in Table 2.

Table 2

The Butcher tableau corresponding to the four stages method.

0
c_2	a21
c_3	a31	a32
c_4	a41	a42	a43
	b1	b2	b3	b4

The well–known RKAM(4,4) technique can be presented in the form of Butcher array [22]. The Butcher tableau, corresponding to this method is presented in Table 3.

Table 3

The Butcher tableau corresponding to the RKAM(4,4) method.

0
12	12
12	0	12
1	0	0	1	0
1	16	13	13	16

Table 3 yields the equivalent corresponding equations defining the RKAM(4,4) method:whereWith Butcher array, the RKAM(4,4) method may also be presented in the modified form shown in Table 4.

Table 4

The Butcher tableau corresponding to the modified form of RKAM(4,4) method.

0
12	0	12
1	0	0	1
	13	13	13

Again, the RKRMS(4,4) method as in [15] can be set to the following form:where

The combination of the schemes, RKAM(4,4) and RKRMS(4,4) given by Eqs. (14) and (15)) is referred to as RKARMS(4,4) scheme and can be expressed as [22].

As in [22], the RKARMS(4,4) algorithm can be presented in the form of Butcher array. The Butcher tableau, corresponding to this method is presented in Table 5.

Table 5

The Butcher tableau corresponding to the RKARMS(4,4) method.

0
12	12
12
	0	12
1	0	0	1
	…	…	…
	…	…	…
12	116	716
1	18	−1756	3328
	13	13	13
	13	13	13
	E^T

Hence by Butcher array (Table 1), one can setand the estimate of the LTE as

In our employed embedded technique, four stages mentioned above are required to attain the approximate solution, which share the set of vectors and using whereas .

Analysis of EC in RKARMS(4,4) method

The basic categories of errors in numerical solutions of ODEs and PDEs are rounding and truncation errors. Here, the truncation error for controlling step size is only discussed.

Derivation of local truncation error for RKARMS(4,4) algorithm

As in Lotkin [28] and Ralston [29], an ERREST for the RKAM(4,4) scheme can be given by , where and are constants such that and . To regulate , Eq. (20) may be used to attain an estimation of the LTE for the embedded RKARMS(4,4) scheme as .

The LTE for the well-established RKAM(4,4) scheme is given by and that for the RKRMS(4,4) is , where and are the approximated field values at attained through the methods RKAM(4,4) and RKRMS(4,4), respectively, whereas and are the respective attained by the RKAM(4,4) and RKRMS(4,4) methods. Therefore, an ERREST for the approximations of the field values at can be given by means of

As in [23], the LTE of the RKAM(4,4) scheme can be given bywhere the LTE of the RKRMS(4,4) scheme is given by

The absolute difference between and can then be set as

As in Lotkin [28], with the intention of controlling the error, the selection of the step size can be made, which leads to [24]where and are constants.

Error estimation for the RKARMS(4,4) method

Here, for the RKARMS(4,4) technique with EC program, the ERREST can be set as the difference between the numerical estimates of the field values attained by the methods that were used to embed. The ERREST then from Eq. (24) can be put to the form:

However, a comparison of the LTE, global truncation error (GTE) and ERREST for the RKARMS(4,4) algorithm with those for the RK–embedded algorithms, namely RKAHeM(4,4) (combination of RKAM(4,4) and RKHeM(4,4) methods) and RKAHM(4,4) (combination of RKAM(4,4) and RKHM(4,4) methods) is shown in Table 6 for better understanding [30]. It can be observed from Table 6 that the ERREST of the embedded method implemented in the study is less over that of the RKAHeM(4,4) and RKAHM(4,4) methods.

Table 6

Comparison of LTE, GTE, ERREST for different fourth order four stages explicit RK–embedded algorithms.

RK–embedded algorithm	Local truncation error (LTE)	Global truncation error (GTE)	Error estimation (ERREST)
RK-embedded root mean square (RKARMS(4,4))	LTE_AM−LTE_RMS≤(6457184320)×P^4.Qh^5=|yn+1AM−yn+1RMS|≤6457184320P^4.Qh^5	|ε_n|≤(h^4184320LD)×M(e^DL(x_n−x_0)−1)	ERREST=|y_AM−y_RMS|×6457184320
RK-embedded Heronian mean (RKAHeM(4,4))	LTE_AM−LTE_HeM≤(1218091658880)×P^4.Qh^5=|yn+1AM−yn+1HeM|≤1218091658880P^4.Qh^5	|ε_n|≤(h^41658880LD)×M(e^DL(x_n−x_0)−1)	ERREST=|y_AM−y_HeM|×1218091658880
RK-embedded Harmonic mean (RKAHM(4,4))	LTE_AM−LTE_HM≤(546969120)×P^4.Qh^5=|yn+1AM−yn+1HM|≤645769120P^4.Qh^5	|ε_n|≤(h^469120LD)×M(e^DL(x_n−x_0)−1)	ERREST=|y_AM−y_HM|×546969120

Results and discussion

Following the code of stellar evolution, we have estimated the initial profiles of some protoplanets formed by gravitational instability confining the range of mass to 0.3–10 M by an embedded RKARMS(4,4) technique under approximated zero BCs. The results came out through our calculations are presented numerically as well as graphically for discussion and testing validation, as well as testing model efficiency.

Figures 1 and 2 depict, respectively, the initial temperature and pressure profiles inside the protoplanets having masses 0.3, 1, 5 and 10 M. It can be perceptible from Figure 1 that massive protoplanets have hotter interiors. On the other hand, Figure 2 shows that as massive as is a protoplanet, so higher is its surface and central pressures. Our estimated distribution of temperature is found to be comparable well with those obtained in [10, 12], and the pressure profiles are also found to be in an agreement with the corresponding findings presented in [12], but our model presents protoplanets having lesser central pressure over those predicted in [9]. Figure 3 illustrates our model simulated initial configurations for density of the protoplanets having assumed masses. It can be perceived from Figure 3 that a massive protoplanet is denser than those with lower masses except the protoplanets having masses 1 and 10, which is found to be consistent with the results obtained in [10, 12]. It is found that the protoplanet having mass 1 is less centrally condensed in comparison with the protoplanet with mass 0.3. In contrast, the protoplanet with mass 10 is found to be less dense in comparison with the protoplanets having mass 5 as well as that of 7, which agrees fairly well with the corresponding findings presented in [12]. In reality, different numerical models have been found to predict different initial profiles inside the protoplanets formed via gravitational instability. For instance, simulations conducted by Boss [8, 31] predict less dense and colder protoplanets than the simulation made by Mayer et al. [32, 33] and the estimations made by the said investigators present denser and warmer initial profiles of the protoplanets than their initial profiles presented by DeCampli and Cameron [9] and Bodenheimer et al. [26]. Nonetheless, our numerical simulation possesses unique solution suggesting that the disk instability mechanism is a reasonable hypothesis.

Figure 1

Temperature profiles inside some initial protoplanets.

Figure 2

Pressure profiles inside some initial protoplanets.

Figure 3

Density profiles inside some initial protoplanets.

Likewise, we have used the RKAM(4,4) and RKRMS(4,4) schemes in estimating the results for making a comparison. The results attained by the RKARMS(4,4) method were found to be comparable with those attained through both the schemes, RKAM(4,4) and RKRMS(4,4). To save space consumption, only the initial profiles of physical variables for a protoplanet having 5 are presented in tabular form in Table 7 for better perspective. From Table 7, it is inferred that the outcomes obtained by the RKARMS(4,4) technique agree fairly well with those attained through the RKAM(4,4) and RKRMS(4,4) techniques, but a change in the outcomes can be found nearer to the core region.

Table 7

Comparison of computed results for varying x (0.99–0.001) by the RKAM(4,4), RKRMS(4,4) and novel RKARMS(4,4) algorithms in the case of estimating thermodynamic variables inside a 5 Jupiter mass protoplanet.

r/R	Classical 4th order RK method			4th order RKRMS method			Novel RKARMS(4,4) algorithm
	P (dyne cm−2)	T(K)	ρ (gm cm−3)	P (dyne cm−2)	T(K)	ρ (gm cm−3)	P (dyne cm−2)	T(K)	ρ (gm cm−3)
0.99	0.0045	8.1253	6.0044e−12	0.0045	8.1253	6.0040e−12	0.0045	8.1253	6.0044e−12
0.9	1.7970	89.0003	2.3842e−10	1.7973	89.0004	2.3832e−10	1.7976	89.0004	2.3837e−10
0.8	12.9651	196.1498	8.5961e−10	12.9602	196.1518	8.5928e−10	12.9626	196.1508	8.5944e−10
0.7	44.4232	321.0118	1.9636e−09	44.4083	321.0213	1.9629e−09	44.4156	321.0167	1.9632e−09
0.6	109.3609	460.2797	3.6254e−09	109.3335	460.3089	3.6245e−09	109.3468	460.2948	3.6249e−09
0.5	218.8593	607.4954	5.8043e−09	218.8342	607.5670	5.8036e−09	218.8464	607.5323	5.8040e−09
0.4	374.5993	753.1886	8.2789e−09	374.6351	753.3408	8.2796e−09	374.6178	753.2670	8.2793e−09
0.3	559.7960	884.4817	1.0620e−08	560.0344	884.7774	1.0624e−08	561.8516	885.8542	1.0643e−08
0.2	747.2811	992.8201	1.2387e−08	748.0520	993.3925	1.2399e−08	747.6783	993.1149	1.2393e−08
0.1	889.9564	1064.6931	1.3112e−08	892.2707	1065.9745	1.3146e−08	891.1485	1065.3533	1.3130e−08
0.01	1048.5194	1136.8615	1.4044e−08	1076.9385	1149.2760	1.4425e−08	1063.0931	1143.2524	1.4239e−08
0.001	2541.6040	1620.7022	3.3736e−08	3123.2788	1760.2253	4.1457e−08	2788.1316	1681.2756	3.7009e−08

With the aim of comparing computational effectiveness of the proposed technique with those of RKAM(4.4) and RKRMS(4,4), the codes produced by the techniques were implemented on the computer having configuration ‘Intel(R) Core(TM) i5-4570, 4th generation’ with different initial time steps and different starting points (0.05, 0.01 and 0.001). The total computational time in each case for the RKARMS(4,4) method was attained to be a little bit more over the methods that were used to generate the embedded RKARMS(4,4) method. In terms of computational cost, the results of our calculation and that of Paul and Senthilkumar [25] are striking. The reason behind may be the choice of the tolerance. Because in the case of the embedded method like RKARMS(4,4), the adoptive step size is advanced in each step of the solution procedure, if necessary, depending upon the accuracy of results. However, in each time step, repeated calculations may be needed until the desired order of accuracy is attained. But in both the RKAM(4,4) and RKRMS(4,4) schemes, the same step size is maintained throughout a given domain of integration. Table 8 clarifies the fact mentioned above that RKARMS(4,4) needs less number of time steps but a little bit more computational time (a fraction of second) over the methods, RKAM(4,4) and RKRMS(4,4). But, the computational cost achieved in the case of the newly proposed method is not a big issue. However, at this juncture, it is justified to note that in a system of ODEs, the unknowns are dependent to each other. Thus, if in the case of an unknown, an error is retained, then it influences the calculations of the unknowns for the next subsequent steps and errors can be piled up as well. In the case of the RKARMS(4,4) method, a result with a certain accuracy can be obtained depending upon setting up the tolerance. But there is not such a type of advantage available in RKAM(4,4) and RKRMS(4,4) techniques. The outcomes achieved by both the RKAM(4,4) and RKRMS(4,4) techniques solely depend on the selection of step size, where there is no way of advancing the step size. This leads to increase the number of total time step which in turn may help increasing computational error. Thus, with regard to accuracy, the method adopted here in this study can be found to defeat the other two methods.

Table 8

Comparison of our computed results for central values calculated by the RKAM(4,4), RKRMS(4,4) and RKARMS(4,4) methods. The calculations here are made for a protoplanet with 5M with different initial time steps. Starting values are different but the calculations are made downwards to the point 0.99 in each case.

Method	Initial time step and starting value	Total step needed	Computational time (second)	P (dyne cm−2) at the end point	T (°K) at the end point
RKAM(4,4)	0.05	9400	0.13	977.41	1105.36
	0.01	9800	0.12	1046.35	1135.91
	0.001	9890	0.11	1048.96	1137.05
RKRMS(4,4)	0.05	9400	0.16	977.86	1105.56
	0.01	9800	0.15	1047.86	1136.56
	0.001	9890	0.15	1073.66	1147.84
RKARMS(4,4)	0.05	3232	0.32	978.70	1105.94
	0.01	3622	0.33	1060.18	1141.98
	0.001	4421	0.41	1635.46	1362.13

For testing efficiency of the method adopted in the study, ERRESTs are calculated with reference to and while calculating them from 0.01 downwards to the point 0.999 in the case of a and are depicted in Figure 4. The figure clearly shows the effective applicability of the different approach adopted in this study. It is to be noted here that there is no opportunity to estimate errors in the case of both the RKAM(4,4) and RKRMS(4,4) techniques.

Figure 4

ERREST in P and T. Here 1M is considered.

Further, for testing the suitability of the methods in solving end point constraints, we have estimated our outcomes with variable end points by the methods taken into account. Significant changes in results at the assumed end points are observed when calculating by the techniques, RKAM(4,4) and RKRMS(4,4) but in the case of the novel embedded RKARMS(4,4) technique, the outcomes are seen to be consistent for the variable end points. The results in the case of a are depicted in Figure 5 for better understanding. Other similar such depiction of figures is not made for the intent of conciseness, which, we think, is also not physically instructive. But the same analyses can be found to make for the other protoplanets with the assumed masses. Therefore, the embedded RKARMS(4,4) technique is obtained as more suitable for solving protoplanetary structure equations over the established RKAM(4,4) method as well as RKRMS(4,4) method with respect to accuracy, effectiveness and solving end point constraints. One shortcoming of this approach is that it may increase the computational overhead significantly depending on the selection of the opening step size as well as the tolerance. But computational cost as obtained by the study cannot be a considerable issue in the present day computers with high performance.

Figure 5

Calculations for P and T with varying end points.

Conclusion and future perspective

Implementation of the anew proposed RKARMS(4,4) algorithm to effectively estimate and validate the initial configurations of temperature, pressure and density inside some protoplanets in their initial stage formed via disk instability for protoplanetary masses between 0.3 and 10 is the main focus of this paper. It is observed from the theoretical point of view that the newly developed RKARMS(4,4) technique has less LTE and less ERREST over the embedded RKAHeM(4,4) and RKAHM(4,4) methods and is found to be more suitable over the RKAM(4,4) and RKRMS(4,4) methods in solving end point constraints. Therefore, the RKARMS(4,4) method can be an alternative efficient technique to study and analyze the evolution of extra-solar giant planets formed by gravitational instability and hence to study and analyze other nonlinear problems arising in engineering, mathematics, physics and other branches of science. Our impending investigation will be concentrated towards analyzing, investigating and employing sequential higher order numerical integration and parallel algorithms in addition to space and time complexity.

Declarations

Author contribution statement

Gour Chandra Paul: Conceived and designed the analysis; Analyzed and interpreted the data; Wrote the paper.

Sukumar Senthilkumar: Analyzed and interpreted the data.

Hafijur Rahman: Analyzed and interpreted the data; Wrote the paper.

Funding statement

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Competing interest statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.