ABSTRACT

In search of nonlinear oscillations, we envision a 3D elliptic curvature-dependent nonuniform charge distribution to creating an electric field along the symmetry axis causing a massive point-like charged particle placed on the symmetry axis to oscillate in a delayed/hesitant nonlinear mode. The charge distribution is a 3D twisted line creating nontrivial electric field causing an unexpected oscillation that is non-orthodox defying the common sense. Calculation of this research flavored investigation is entirely based on utilities accompanied with Computer Algebra Systems (CAS) especially Mathematic [1]. The characteristics of the delayed oscillations in addition to embodying classic graphics displaying the time-dependent kinematic quantities are augmented including various phase diagrams signifying the nonlinear oscillations. The output of our investigation is compared to nonlinear non-delayed oscillations revealing fresh insight. For comprehensive understanding of the hesitant oscillator a simulation program is crafted clarifying visually the scenario on hand.

Keywords:

Nonlinear Oscillator, 3D Elliptic Curvature-Dependent Charge Distribution, Computer Algebra System, Mathematica

1. Introduction, Motivations and Goals

In our previous work [2] we showed how the eccentricities of elliptic curvature-dependent charge distribution impact the characteristics of nonlinear oscillations of a massive point-like charged particle. The main physical motivation was to handle the nonuniformity of the charge distribution resulting from the elliptical curvature. As explained, we established the conceptual similarities between the notion of the mechanical stress of an ellipse and the elliptic charge distribution. In short, for a 2D ellipse the smaller the curvature radius the sharper the bent and hence higher the mechanical stress. And from charge distribution point of view, the sharper the bent the higher the charge density and vice versa. Calculation of the needed electric field created by such nonuniform charge distribution encounters challenging integration with hopeless analytic long-hand and failed Computer Algebra Systems [1] [3]. Utilizing the numeric features of CAS, we overcome the challenges introducing a semi-numeric-analytic computational method. This method not only solves the issue on hand it has the flexibilities to be applicable to similar scenarios. Reference [2] detailed descriptive explanation of our approach accompanied with robust Mathematica codes.

Our current investigation is a derivative, an upgrade of our previous work. It is upgraded two-folds. We consider 1) a 3D elliptic contour housing the charge distribution and, reference [2] dealt with 2D untwisted distribution. 2) the contour is wrinkled including several ripples. The curvature of the 2D ripple-free ellipse is given by a continuous varying curvature expressed by a relatively simple cyclic function as shown in Figure 1, [2]. Major mathematical and computational challenges all are addressed in [2], interested readers are invited to review [2]. For the current scenario, the curvature function is not trivial. More on this in follow up paragraphs.

With these upgrades it is the goal of our research flavored investigation to investigate the impact of the 3D elliptic contoured frame housing the charge distribution on the motion of a massive point-like charged particle placed on the vertical symmetry axis of the contour. It is intuitive to foresee the impact. Having the charge on the contour to have the same “sign” as the charge of the point charge the mutual electrostatic repulsive force between them may be adjusted to compete against the gravity pull. These two oppositely oriented forces for suitably chosen physical quantities such as mass, charge, size of the contour etc. may result oscillations. The nature of the electrostatic potential and its related electric field should qualify nonregular, nonlinear oscillations. Calculations and crafted computer codes quantitatively justify our hinge.

With this goal we craft this report, it is composed of three sections. In addition to Section 1, Introduction, Motivations and Goals in Section 2 we explain the physics of the problem laying-down the mathematical formulation. Since the objective is to apply CAS, we include the needed Mathematica codes. Utilizing the output, we then display the relevant kinematic time-dependent quantities. As pointed out in the Abstract in addition to the classic graphs we graph various phase diagrams assisting to comprehension of the nature of the nonlinear, so suitably named delayed/hesitant nonlinear oscillations. The last section is the Discussions and Conclusions.

2. Physics of the Problem and Formulation

We begin with displaying the scenario on hand. Figure 2 shows a segment of an elliptic-based cylinder. The contour drawn on the shell is an elliptic curve. Meaning from the top view looking down the symmetry axis the contour likes looks an ellipse. The code generated the shown snapshot has the capability of changing the view making the explained feature visual. From the side view the contour would appear as a pair of six humped-valley curve validating visually the 3D feature of the curve. The contour houses the charge. In the shown figure the number of the humps and valleys as an example is set twelve; in practice this can easily be set to meet the need. The dot on the symmetry axis is a typical point of interest where a loose point-like charge is eventually placed. The elliptic curvature of the contour distributes the charge accordingly making the distribution cyclically uneven. To calculate the electrostatic force that the contoured charge distribution exerts on the dot there is a need to calculate the electrostatic

field. As we had somewhat encountered a similar situation [2] first we calculate the electrostatic potential and then derive its associated electric field. In doing so nonetheless, we encounter an integration of an integrand that depends on the inverse length of the shown vector. The tail of the vector slides along the rim of the contour and its head points to the shown dot on the symmetry axis. By integrating the inverse length of the arrow, the integration sweeps and collects the contributions of the differential charges. This feature is built in the shown animation where the slider φ runs over one complete cycle. By clicking the (+) the code runs automatically displaying the revolving arrow. The animation shows the changing φ impacts the length of the 3D vector and that adds to the challenges evaluating the integral.

Calculations and Codes

Because the charge distribution on the 3D contour shown in Figure 2 is curvature-dependent, it is not uniform, first we identify the curvature function for later use. This function is given by [4],

$\text{curvature}=\frac{\left|r^{\prime }\times r^{″}\right|}{\left|{r^{\prime }}^3\right|}$

where $r$ is the 3D vector with its tail at the origin and its head touching the contour. The prim and double-prime are the first and the second derivates of the position vector, $r$ with respect the variable parameter within the $r$, respectively. Vector $r$ is chosen as, $r=\left\{a\cos \left[\phi \right],b\sin \left[\phi \right],c\cos \left[n\phi \right]\right\}$. Here, a and b are the semi major and minor axes lengths of the ellipse on the horizontal x-y plane. The c and n are the vertical amplitude of the contour along the z-axis, and n controls the number of ripples, respectively. Accordingly, the shown contour in Figure 2 corresponds to $\left\{a,b,c,n\right\}=\left\{2.,0.5,0.5,6\right\}$. The eccentricity of the ellipse is 0.96 representing a pinched ellipse along the y-axis. Such eccentricity makes the ellipse to deviate being looks like a circle, pronouncing its physical impact. The azimuthal angle is, $0\le \phi \le 2\pi$ constituting the explicit variable

parameter allowing to trace the contour. Hence, $r^{\prime }=\frac{\text{d}}{\text{d}\phi }r\left(\phi \right)$ and $r^{″}=\frac{{\text{d}}^2}{\text{d}{\phi }^2}r\left(\phi \right)$. The code needed to calculate the curvature function is,

r1[a_, b_, c_, n_, φ_] = {a Cos[φ], b Sin[φ], c Cos[n φ]};

rp[φ] = D[r1[a, b, c, n, φ], {φ, 1}];

rpp[φ] = D[r1[a, b, c, n, φ], {φ, 2}];

curvature[a_, b_, c_, n_, φ] = sqrt((rp[φ] × rpp[φ]).(rp[φ] × rpp[φ])) /(sqrt((rp[φ]).(rp[φ])))³//Simplify;

Plot[Evaluate[curvature[2, 0.5, 0.5, 6, φ]], {φ, 0, 2π}, Ticks → {Range[0, 2π, π/3], All}, GridLines → {Range[0, 2π, π/3], Automatic}, AxesLabel → {“φ”, “Curvature”}, Plot Range → All, PlotStyle → Black], (1)

As shown the curvature function and hence the charge distribution on the contour is nonuniform. The spikes are indicatives of the sharp twists, i.e., the sharper the curvature the high the concentrated charges, contrary to the valleys and the segments within. Noting there is a substantial difference between the character of the curvature function in 3D and 2D [2]. The twists and turns of the 3D contour are embodied within the curvature function; pictorially justifying that the electric field along the contour’s vertical symmetry axis would differ from its 2D flat elliptic case [2]. As one anticipates the ripples on the rim of the 3D contour, its curvature function and associated charge distribution would affect the electric field and hence the electric force along the symmetry axis resulting potential oscillations other than harmonic. This has been quantitatively confirmed in the following paragraphs.

As we experienced [elliptic charge] the “straight forward” traditional calculation of the electric field for a class of unorthodox charge distributions is unsuccessful. Computer Algebra Systems such as, Mathematica [1] and Maple [5] are incapable of producing useful output; the current investigation falls in the same category. Henceforth, we follow the same successful method we introduced [ellipse]. First with the ultimate goal of calculating the electric field and hence its associated force, $F=qE$, with $E=-\nabla V$ we evaluate the electrostatic potential, V numerically and then by fitting the tabulated data generate a continuous function evaluating the $E$. Here the potential V is,

$V\left(z\right)={\displaystyle {\oint }_C\lambda \left(\phi \right)\frac{1}{\text{distance}\left(z,\phi \right)}\text{d}\mathcal{l}}$ (2)

The first term of the integrand in (2) is the linear-charge density that explicitly depends to the azimuthal angle φ curling about the contour, the second term is the symbolic generic function indicating the inverse distance character of the potential. The “distance(z, φ)” is the distance from the differential charge element on the rim of the contour to a point on the symmetry, z-axis. The integration is performed on the contour, C.

Here are the quantities of interest. The length of the contour, this is evaluated numerically,

curve length = N Integrate[sqrt(rp[φ].rp[φ])/.{a → 2, b → 0.5, c → 0.5, n → 6}, {φ, 0, 2π}] (3)

15.47

As expected, the contour length of the rippled 3D ellipse, 15.47 is longer than the circumference of the 2D ellipse with the same geometric characters. The perimeter of the unrippled corresponding ellipse is calculated directly via numeric integration on one hand and via applying Ramanujan’s [6] approximation on the other hand; these values are the same and are,

{elliptic Length, Ramanja} = {N Integrate[sqrt(rp[φ].rp[φ])/.{a → 2, b → 0.5, c → 0., n → 6}, {φ, 0, 2π}], π(3(a + b)-sqrt((3a + b)(a + 3b)))/.{a → 2, b → 0.5}} (4)

{8.57, 8.57}

We then calculate the $\text{distance}\left(z,\phi \right)$. Two vectors are needed, $r_1$ and $r_2$ and the distance is $\left|r_1-r_2\right|$. These are given via codes,

r2[z_] = {0, 0, z}

r1[a, b, c, n, φ]-r2[z]

dist[φ_, z_] = sqrt((r1[a, b, c, n, φ]-r2[z]).(r1[a, b, c, n, φ]-r2[z])) (5)

Noting the curvature function has a dimension of 1/L. To make it unitless (dimensionless) we multiply by 1/2(a + b) i.e., half the average size of the ellipse.

Since the linear-charge density $\lambda \left(\phi \right)$ ought to be normalized i.e., ${\displaystyle {\oint }_C\lambda \left(\phi \right)\text{d}\mathcal{l}}=q_1$, where q₁ is the charge on the contour, we introduce a normalization factor, it is given via the code,

normalization Factor = N Integrate[1/curve length 1/2 (a + b)curvature[a, b, c, n, φ] sqrt(rp[φ].rp[φ])/.{a → 2, b → 0.5, c → 0.5, n → 6},{φ, 0, 2π}], (6)

2.279

To pursue the needed calculations, we define the needed practical physical parameters via, a list, values. More explanation of m and g in the forthcoming paragraphs.

values = {k → 9 × 10^9, q1 → 1 × 10^−6, q2 → 2 × 10^−6, m → 0.1 × 10^−3, g → 9.8} (7)

The units are SI, the $k=\frac{1}{4\pi {ϵ}_0}=9\times {10}^9$, with ϵ0 being the permittivity of vacuum. The q₁ and q₂ are the charges on the contour and the secondary point-like particle, respectively. The m and the g are the mass of the particle and the gravity acceleration, respectively. Putting all this together the code for the charge density and the integrand of Equation (1) reads,

λ[φ_] = 1/normalization Factor q1/curve length 1/2(a + b)curvature[a, b, c, n, φ]/.values potential integrand[z_, φ_] = λ[φ]1/(dist[φ, z]) sqrt(rp[φ].rp[φ])/.{a → 2, b → 0.5, c → 0.5, n → 6}//Simplify; (8)

Next, we calculate the electrostatic potential V(z) for a set of z-values. This is accomplished via numeric integration. The output is displayed, Figure 3.

V[z_] = Table[{z, N Integrate[Evaluate[k potential integrand[z,φ]/.{a → 2, b → 0.5, c → 0.5, n → 6}]/.values, {φ, 0, 2π}]}, {z, 0, 10, 1}];

listplot V[z] = ListPlot[V[z], AxesLabel → {“z(m)”, “V(z)Volt”}, GridLines → Automatic, PlotStyle → Black] (9)

The shown data appears somewhat exponential. To calculate the needed electric field, $E$ we search for a one variable analytic function by fitting the data shown in Figure 3. The steps, graphs and compared fit to data are coded, see Figure 4.

modelz[z_]: = c1 + d1 z + e1 e^(f1z^2)

fitVz = FindFit[V[z], modelz[z], {c1, d1, (*g1,*)e1, f1}, z];

plotfitdataz = Plot[modelz[z]/.fitVz, {z, 0, 10}, PlotStyle → Black, AxesOrigin → {0, 0}];

s = Show[listplotV[z], plotfitdataz] (10)

By more aggressive search one may find a “refined” function fitting the data; for time being we are satisfied with the fit. The proposed model function is compatible with our reference work [2]. With this function in hand, we pursue calculating the electric field and the electrostatic force exerting on the loose massive point-like charge, q₂.

Efield = -D[Evaluate[modelz[z]/.fitVz],{z,1}];

Eforce = q2 Efield/.values//Simplify;

Plot[10^3Eforce/.values, {z, 0, 10}, AxesOrigin → {0, 0}, PlotRange → All, AxesLabel → {“z(m)”, “Eforce(mN)”}, GridLines → Automatic, PlotStyle → Black] (11)

Noticing the force along the symmetry axis has a humped shape within the length of the ellipse’s major axis, for higher heights it stays constant. One therefore anticipates that the interesting physics phenomena should occur if a loose point-like charge is placed within the humped region. This has been justified furthering calculations. The practical units of the force shown in Figure 5 is mN guided with this value and by trial and error we search for a suitable mass, this is embedded in the previously defined list, values. The ordinate of the equilibrium where the electrostatic force compensates the gravity force comes about solving $F_{net}=0$.

Solve[Evaluate[Eforce/.values]==mg/.values, z], (12)

{{z → 0.132}, {z → 3.43}}

These two values fall within the desired height. Guided by these values, we place the q₂ within the vicinity of one of these points, e.g., the lower height. This sets up one of the initial conditions to solving the equation of motion. Applying Newton’s law, $F_{net}=ma$ the equation of motion, its numeric solution and graphic output are,

EquationOfMotion = z”[t]-(Eforce/.z → z[t])((1/m)/.values) + g/.values

solEquationOfMotion = NDSolve[{EquationOfMotion==0, z[0]==0.135, z'[0]==0}, z[t], {t, 0, 20}]

plotingOscillations = Plot[Evaluate[z[t]/.solEquationOfMotion], {t, 0, 20}, AxesLabel → {“t(s)”, “z(m)”}, PlotStyle → Black, GridLines → Automatic, AxesOrigin → {0, 0}, PlotRange → All]; (13)

This graph shows the impact of the charge distribution on the uncommon oscillations of the free point-like charge. It is a vivid signature of nonharmonic oscillations. Figure 6 shows according to applied initial conditions, q₂ is released freely from the vicinity of the lower equilibrium position along the z-axis. At the beginning it is pushed away from the contour, reaches the maximum height, it then after a “split” second it falls back to its initial position. At the bottom of its fall contrary to the maximum height scenario it appears it “takes it time” relaxing for a good one second and then builds up energy with “hesitation” jumping back to the max height; in the absence of friction, it repeats the process all over again. To make the scenario comprehensible we include velocity and acceleration of the particle vs. time as well. These are collectively shown in Figure 7.

This figure is descriptive. As such the velocity shown in short dashed Gray describes the particle starting from rest because of the electrostatic force gains speed, after reaching to its max speed on its way to the max height it slows. At the max height it stops for a “split” second reverses direction heads to the initial position. At the bottom it sits there for a while, as if it is motionless, builds up energy jumping back up, repeating the oscillations. Acceleration shown in long dashed Gray is similarly descriptive, it is left to be justified by the reader.

For the sake of completeness additionally we include two phase diagrams, they help to the understanding of the impact of the nonuniform charge distribution on the characters of the “hesitant” oscillator, see Figure 8.

The interested readers are encouraged to compare the specifics of the shown diagrams to the phase diagrams of harmonic, linear oscillations, they are drastically different.

3. Discussions and Conclusions

Motivated with our latest work [2] it was naturally essential to augment our investigation considering a 3D elliptic curvature-dependent charge distribution anticipating conducive nonlinear oscillations. To quantifying the mathematical challenges such as calculation of the electrostatic force we justified that the routine classic approach was fruitless. We showed also direct application of the computational utilities of the Computer Algebra Systems (CAS) have short-comings and CPU expensive. Henceforth, we applied our semi-numeric-symbolic method circumventing the mathematical challenges conducive to fruitful calculation of the needed physical quantities. Equation of motion resulting from the force is nonlinear differential equation with no analytic solution. We solve the equation applying Mathematica’s numeric routine. Having these solutions on hand by displaying the kinematic time-dependent quantities such as position and velocity vs. time we showed the meaning of the unorthodox nonlinear oscillations, giving a meaning to the suitably chosen title “hesitant” oscillator. The last section of our work embodies suitably crafted phase diagrams assisting to the understanding of the delayed nonlinear oscillations.

Interested readers may utilize formulation and codes embedded in our work to augmenting the scope of the investigation. Projects such as the impact of elliptic eccentricities, variation of the number of the elliptic ripples, and initial conditions on the characteristics of the oscillations have the potential of furthering the scope of our work. Readers may also find [1] [7] to creating graphs resourceful.

Acknowledgements

The author gracefully acknowledges the John T. and Paige S. Smith Professorship funds for completing and publishing this work.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

Cite this paper

Sarafian, H. (2021) Nonlinear Electrostatic “Hesitant” Oscillator. American Journal of Computational Mathematics, 11, 42-52. https://doi.org/10.4236/ajcm.2021.111004

References