The modified homotopy perturbation method with an auxiliary term for the nonlinear oscillator with discontinuity

Abstract

Recently, a modified homotopy perturbation method with an auxiliary term was proposed for obtaining an analytical approximate solution of a nonlinear equation. The auxiliary term involves a tuning parameter that gives the modified homotopy perturbation method an additional degree-of-freedom to construct the approximate solution. In this paper, we present a systematic approach for choosing the auxiliary parameter and furthermore show that an optimal value of the parameter yields significant improvement in the approximate solution compared to the standard homotopy perturbation method. We consider the nonlinear oscillator with a discontinuity to illustrate the proposed approach.

Introduction

Nonlinear differential equations play a fundamental role in physics, applied mathematics, engineering, and other disciplines. The last two decades have seen a surge of interest in analytical tools for studying such equations. In particular, the nonlinear oscillator has received considerable attention where the aim is to jointly approximate the frequency of oscillation and the solution of a nonlinear equation. This problem is further compounded by the presence of a discontinuous term in the oscillator equation.

The absence of a small parameter precludes application of the classic perturbation methods. ¹ We begin a survey ^a by reviewing a class of methods that leads to a rapid solution with reasonably good accuracy. This class of methods includes, for example, the amplitude–frequency formulation^2–4 and the max–min approach. ⁵ The application of the amplitude–frequency formulation to a nonlinear oscillator with a discontinuity was discussed by Zhang. ⁶ Zeng et al. ⁷ carried out an empirical study of the amplitude–frequency relationship for the same problem. Recent developments to the approach can be found in the work of He.^8,9 Zeng ¹⁰ studied the amplitude–frequency relationship using the max–min approach for a nonlinear oscillator with a discontinuity. The periodic solution of a generalized nonlinear discontinuity equation using the max–min approach was discussed by Demirbağ and Kaya. ¹¹

Next, we review a class of iterative methods where higher order approximations can deliver more accurate solutions. The approaches in this class include, for example, the homotopy perturbation method (HPM),^12–15 the modified Lindstedt–Poincaré method,^16–20 the parameter-expansion method,^21–24 energy balance method,^25,26 and the variational iteration method.^12,27,28 The application of the HPM to a nonlinear oscillator equation with a discontinuity was discussed for the first time by He. ¹⁴ Özis and Yildirim ²⁹ showed that the construction of the homotopy equation is problem dependent and affects the accuracy of approximation. In a series of papers,^17,18,30 He showed that modified Lindstedt–Poincaré methods could be applied to strongly nonlinear equations. Liu ¹⁹ applied He’s method ¹⁷ to approximate the period of a nonlinear oscillator with a discontinuity. Wang and He ²³ and Zengin et al. ²⁴ demonstrated the effectiveness of the parameter-expansion methods for nonlinear oscillators with discontinuities. Application of the energy balance method to nonlinear oscillators with discontinuities was discussed by Afrouzi et al. ³¹ and Zhang et al. ³² The variational iteration method for a nonlinear oscillator was discussed by He. ³³ Rafei et al. ³⁴ applied the approach to a nonlinear oscillator with a discontinuity.

In all of the above-mentioned works only first-order approximate solutions were constructed which were found to be in good agreement with those obtained via competing methods. Beléndez et al.^35–37 considered an anti-symmetric quadratic nonlinear oscillator, ³⁸ a nonlinear oscillator with a discontinuity, and developed higher-order analytical approximate solutions using a modified HPM and the generalized harmonic balance method. Kaya and Demirbağ ³⁹ discussed a class of nonlinear discontinuity equations which subsumes the nonlinear oscillators considered in the above-mentioned works. Furthermore, the authors applied the parameter expanding methods to construct first-order approximate solutions that agree well with those obtained via competing methods.

A number of modifications to the standard HPM procedure have extended the applicability of the approach to an even broader class of problems which include the Riccati equation, ⁴⁰ integral equations,^41,42 heat transfer problem, ⁴³ differential-difference equations, ⁴⁴ and optimal control problems. ⁴⁵ Recently, the introduction of an auxiliary term in the homotopy equation has been discussed.^46,47 Noor ⁴⁶ proposed using an auxiliary operator that needs to be determined. Moreover, an auxiliary parameter value of unity was discussed. He ⁴⁷ assumed a particular form of the auxiliary term. For a first-order approximation, the auxiliary term is dropped by setting the auxiliary parameter to zero, reducing the homotopy equation to the standard HPM form. For higher-order approximation, the auxiliary parameter value of unity is expanded as a power series in the homotopy parameter. In both the methods,^46,47 the value of the auxiliary parameter was chosen in an ad hoc manner. In this paper, we discuss for the first time, a systematic procedure for choosing the auxiliary parameter in an optimal manner. We illustrate the proposed approach for the nonlinear discontinuity equation previously considered by Kaya and Demirbağ. ³⁹ As already mentioned, the generalized nonlinear equation encompasses a broad class of strongly nonlinear and nonsmooth problems considered in the literature. We compare the performance of the proposed approach with the standard HPM for obtaining analytic first- and second-order approximation of the frequency of oscillation and the solution of the nonlinear equation.

Modified HPM with an auxiliary term

Consider the general nonlinear differential equation

$$\mathcal{A}(v)-f(t)=0,t\in \Omega$$

(1)

with initial conditions $v(0)=A,v^{\prime }(0)=0$ where $\mathcal{A}$ is a general differential operator and f(t) is a known analytic function. Equation (1) can be written in the form

$$\mathcal{L}(v)+\mathcal{N}(v)-f(t)=0$$

where $\mathcal{L}$ and $\mathcal{N}$ are the linear and nonlinear operators, respectively.

The HPM ⁴⁸ constructs a homotopy of equation (1), $v(t,p):\Omega \times [0,1]\to \mathbb{ℝ}$ satisfying

$$H(v,p)=(1-p)[\mathcal{L}(v)-\mathcal{L}(v_o)]+p[\mathcal{L}(v)+\mathcal{N}(v)-f(t)]=0$$

(2)

where $p\in [0,1]$ is an embedding parameter and v_o is an initial approximation satisfying the initial conditions.

From equation (2), we have

$$H(v,0)=\mathcal{L}(v)-\mathcal{L}(v_o)=0$$

(3)

$$H(v,1)=\mathcal{L}(v)+\mathcal{N}(v)-f(t)=0$$

(4)

As p increases from 0 to 1, the trivial problem (equation (3)) continuously deforms to the original problem (equation (4)).

The modified HPM ⁴⁷ introduced an auxiliary term in the homotopy equation

$$\tilde{H}(v,p)=H(v,p)+kp(1-p)v=0$$

(5)

where k is an auxiliary parameter. For k = 0, equation (5) reduces to the standard homotopy (equation (2)). Note that $\tilde{H}(v,0)=H(v,0)$ and $\tilde{H}(v,1)=H(v,1)$ since the auxiliary term vanishes for $p\in \{0,1\}$ . The additional degree-of-freedom due to the auxiliary parameter k can yield a considerable improvement in the accuracy of the solution series if k is chosen in an optimal manner.

Assuming the solution of equation (1) can be expressed as a power series in p, we have

$$v=v_o+p{v}_1+p^2{v}_2+\mathrm{\cdots }$$

(6)

Furthermore, we assume that the auxiliary parameter k can also be expressed as a power series in p to yield

$$k=k_o+p{k}_1+p^2{k}_2+\mathrm{\cdots }$$

(7)

Substituting equations (6) and (7) in equation (5), we obtain

$$\begin{array}{l}(1-p)[{\displaystyle \sum _{i=0}^{\infty }\mathcal{L}}(p^i{v}_i)-\mathcal{L}(v_o)]+p[{\displaystyle \sum _{i=0}^{\infty }\mathcal{L}}(p^i{v}_i)+\mathcal{N}({\displaystyle \sum _{i=0}^{\infty }{p}^i}{v}_i)-f(t)]+p(1-p)\left({\displaystyle \sum _{i=0}^{\infty }{p}^i}{k}_i\right){\displaystyle \sum _{i=0}^{\infty }{p}^i}{v}_i=0\end{array}$$

Equating coefficients of corresponding powers of p gives a series of linear differential equations with initial conditions. Then, an approximate solution of equation (1) is given by

$$v=\underset{p\to 1}{\lim }v=v_o+v_1+v_2\mathrm{\cdots }$$

The convergence of the series was discussed by He. ¹²

Modified HPM for a discontinuous oscillator

Consider the nonlinear discontinuity equation (6) ³⁹ with $h(v)=\alpha v$

$$v^{\prime \prime }+\alpha v+\lambda \left|v\right|v^{2n-1}=0$$

(8)

for $n\in \mathbb{\mathbb{N}}$ with initial conditions $v(0)=A$ and $v^{\prime }(0)=0$ .

For $\lambda$ not restricted to a small parameter value, classic perturbation methods ¹ are inapplicable. Here, we develop a procedure based on a modified HPM with an auxiliary term ⁴⁷ to construct analytical approximate forms of the frequency of nonlinear oscillation and the solution of equation (8).

Consider the homotopy equation

$$v^{\prime \prime }+\alpha v=-p\lambda \left|v\right|v^{2n-1}+kp(1-p)v$$

(9)

where $p\in [0,1]$ is the homotopy parameter. Note that for p = 0, equation (9) reduces to a linear differential equation

$$v^{\prime \prime }+\alpha v=0$$

which admits an exact periodic solution. For p = 1, the original equation (8) is recovered.

Note that for the auxiliary parameter k = 0, equation (9) reduces to the standard form and for $p\in \{0,1\}$ , the auxiliary term vanishes.

We proceed by expanding the solution v(t), square of the unknown angular frequency $\omega$ and the auxiliary parameter k in equation (9) as a series of the homotopy parameter p

$$v(t)=v_o(t)+p{v}_1(t)+p^2{v}_2(t)+\mathrm{\cdots }$$

(10)

$$\alpha ={\omega }^2-p{\alpha }_1-p^2{\alpha }_2-\mathrm{\cdots }$$

(11)

$$k=k_o+p{k}_1+p^2{k}_2+\mathrm{\cdots }$$

(12)

where ${\alpha }_j,j=1,2,\mathrm{\dots }$ and $k_j,j=0,1,2,\mathrm{\dots }$ will need to be determined.^21,49 For a first-order approximation, He ⁴⁷ takes k_o = 0. Moreover, the identification of $k_j,j=1,2,\mathrm{\dots }$ for higher-order approximation is not discussed. Here, we outline a procedure for choosing the components of k in an optimal manner.

Substituting equations (10) to (12) in equation (9) and equating coefficients of corresponding powers of p on both sides gives a series of linear equations. The first three terms of the series are

$$v_o^{\prime \prime }+{\omega }^2{v}_o=0,v_o(0)=A,v_o^{\prime }(0)=0$$

(13)

$$v_1^{\prime \prime }+{\omega }^2{v}_1={\alpha }_1{v}_o-\lambda \left|v_o\right|v_o^{2n-1}+k_o{v}_o,v_1(0)=0,v_1{(0)}^{\prime }=0$$

(14)

$$v_2^{\prime \prime }+{\omega }^2{v}_2={\alpha }_1{v}_1+{\alpha }_2{v}_o-2n\lambda \left|v_o\right|v_o^{2n-2}{v}_1+k_1{v}_o+k_o{v}_1-k_o{v}_o,v_2(0)=0,v_2^{\prime }(0)=0$$

(15)

where for $h(v)=\left|v\right|v^{2n-1}$ we have used

$$\begin{array}{l}h(v)=h(v_o(t)+p{v}_1(t)+p^2{v}_2(t)+\mathrm{\cdots })\\ =h(v_o)+p{h}^{\prime }(v_o)v_1+p^2[h^{\prime }(v_o)v_2+\frac{1}{2}{h}^{\prime \prime }(v_o)v_1^2]+O(p^3)\end{array}$$

with $h\prime (v)=2n\left|v\right|v^{2n-2}$ and $h″(v)=2n(2n-1)\left|v\right|v^{2n-3}$ .

The solution of equation (13) is

$$v_o(t)=A\cos (\omega t)$$

(16)

Substituting in equation (14), we obtain

$$v_1^{\prime \prime }+{\omega }^2{v}_1=({\alpha }_1+k_o)A\cos (\omega t)-\lambda A^{2n}|\cos (\omega t)|{\cos }^{2n-1}(\omega t)$$

(17)

We take the approach of Beléndez et al.^35,36 and consider a Fourier series expansion of $|\cos (\omega t)|{\cos }^{2n-1}(\omega t)$

$$|\cos \varphi |{\cos }^{2n-1}\varphi ={\displaystyle \sum _{s=0}^{\infty }{b}_{2s+1}}\cos ((2s+1)\varphi )$$

(18)

where

$$b_{2s+1}=\frac{4}{\pi }{\displaystyle {\int }_0^{\pi /2}|}\cos \varphi |{\cos }^{2n-1}(\varphi )\cos ((2s+1)\varphi )\mathrm{d}\varphi$$

(19)

The first term is given by

$$\begin{array}{c}{b}_1=\frac{4}{\pi }{\displaystyle {\int }_0^{\pi /2}|}\cos \varphi |\cos 2n-1(\varphi )\cos (\varphi )\mathrm{d}\varphi =\frac{2}{\sqrt{\pi }}\frac{\Gamma (n+1)}{\Gamma (n+3/2)}\end{array}$$

(20)

Substituting equation (18) in equation (17) we obtain

$$v_1^{\prime \prime }\text{}+{\omega }^2{v}_1\text{}=({\alpha }_{1}A-\lambda b_1{A}^{2n}\text{}+k_{o}A)\cos (\omega t)-\lambda A^{2n}{\displaystyle \sum _{s=1}^{\infty }{b}_{2s+1}}\cos ((2s+1)\omega t)$$

(21)

Dropping the secular terms in v₁ requires

$${\alpha }_{1}A-\lambda b_1{A}^{2n}+k_{o}A=0$$

(22)

Setting p = 1 in equation (11) and substituting for ${\alpha }_1$ from equation (22), the first-order approximate angular frequency is given by

$${\omega }_1=\sqrt{\alpha +\lambda b_1{A}^{2n-1}-k_o}$$

(23)

Note that for k_o = 0 this is the same expression for the first-order approximate angular frequency obtained using the parameter expanding method. ³⁴

With equation (22) satisfied, equation (21) reduces to

$$v_1^{\prime \prime }+{\omega }^2{v}_1=-\lambda A^{2n}{\displaystyle \sum _{s=1}^{\infty }{b}_{2s+1}}\cos ((2s+1)\omega t)$$

(24)

with initial conditions $v_1(0)=0$ and $v_1^{\prime }(0)=0$ . The solution of equation (24) takes the form

$$v_1(t)={\displaystyle \sum _{s=0}^{\infty }{c}_{2s+1}}\cos ((2s+1)\omega t)$$

(25)

Substituting on the left hand side of equation (24)

$${\omega }^2{\displaystyle \sum _{s=0}^{\infty }4}s(s+1)c_{2s+1}\cos ((2s+1)\omega t)=\lambda A^{2n}{\displaystyle \sum _{s=1}^{\infty }{b}_{2s+1}}\cos ((2s+1)\omega t)$$

(26)

Equating coefficients of corresponding harmonics on both sides, the coefficients $c_{2s+1}$ are given by

$$c_{2s+1}=\frac{\lambda A^{2n}{b}_{2s+1}}{4s(s+1){\omega }^2},s=1,2,3,\mathrm{\dots }$$

(27)

The initial condition $v_1(0)=0$ delivers

$$c_1=-{\displaystyle \sum _{s=1}^{\infty }{c}_{2s+1}}$$

(28)

The second-order approximation requires substituting equation (25) in equation (15). But equation (25) has infinite harmonics. Following the approach of Beléndez et al.,^35,36 we truncate the series expansion as

$$v_1^{(M)}=c_1^{(M)}\cos \omega t+{\displaystyle \sum _{s=1}^M{c}_{2s+1}}\cos ((2s+1)\omega t)$$

with $c_1^{(M)}=-{\displaystyle {\sum }_{s=1}^M{c}_{2s+1}}$ .

Note that $\left|c_{2s+1}\right|,s=0,1,2,\mathrm{\dots }$ is monotonically decreasing; only a finite number of harmonics are needed to yield a good approximation.

Consider the case M = 1. Then

$$v_1^{(1)}=c_1^{(1)}\cos \omega t+c_3\cos 3\omega t$$

with $c_1^{(1)}=-c_3$ . Substituting above, we obtain

$$v_1^{(1)}(t)=c_3(\cos 3\omega t-\cos \omega t)$$

(29)

where

$$c_3=\frac{\lambda A^{2n}{b}_3}{8{\omega }^2},b_3=\frac{(n-1/2)}{(n+3/2)}{b}_1$$

(30)

The first-order approximate solution of equation (8) is then given by

$$v_1(t)=v_o(t,{\omega }_1)+v_1^{(1)}(t,{\omega }_1)=(A-c_3)\cos {\omega }_{1}t+c_3\cos 3{\omega }_{1}t$$

(31)

where ${\omega }_1$ is given by equation (23).

Note that ${\omega }_1$ depends nonlinearly on k_o. For fixed values of $n,\lambda ,\alpha$ , and A, an optimal value of k_o can be determined by performing a line search for the minimizer of an appropriate convex criterion such as the squared error in the first-order approximate solution. The numerical solution may be used as a surrogate for the exact solution in case it is not known.

Continuing, we substitute equations (16) and (29) in equation (15) to obtain

$$\begin{array}{l}{v}_2^{\prime \prime }+{\omega }^2{v}_2=({\alpha }_1+k_o)c_3(\cos (3\omega t)-\cos (\omega t))-2n{c}_3\lambda A^{2n-1}|\cos (\omega t)|\cos 2n-2\\ \times (\omega t)(\cos (3\omega t)-\cos (\omega t))+\left({\alpha }_2+k_1-k_o\right)A\cos \left(\omega t\right)\end{array}$$

(32)

Consider the Fourier series expansion

$$|\cos \varphi |{\cos }^{2n-2}\varphi (\cos (3\varphi )-\cos \varphi )={\displaystyle \sum _{s=0}^{\infty }{d}_{2s+1}}\cos ((2s+1)\varphi )$$

(33)

where

$$d_{2s+1}=\frac{4}{\pi }{\displaystyle {\int }_0^{\pi /2}|}\cos \varphi |{\cos }^{2n-2}\varphi (\cos 3\varphi -\cos \varphi )\cos ((2s+1)\varphi )\text{d}\varphi$$

(34)

The first term is given by

$$\begin{array}{l}{d}_1=\frac{4}{\pi }{\displaystyle {\int }_0^{\pi /2}|}\cos \varphi |{\cos }^{2n-2}\varphi (\cos 3\varphi -\cos \varphi )\cos (\varphi )\text{d}\varphi =b_3-b_1\\ =-\frac{2}{n+3/2}{b}_1=-\frac{4}{\sqrt{\pi }(n+3/2)}\frac{\Gamma (n+1)}{\Gamma (n+3/2)}\end{array}$$

(35)

Substituting equation (33) in equation (32) gives

$$\begin{array}{c}{v}_2^{\prime \prime }+{\omega }^2{v}_2=(-{\alpha }_1{c}_3+{\alpha }_{2}A-2n{c}_3{d}_1\lambda A^{2n-1}-k_o{c}_3-k_{o}A\text{}+k_{1}A)\cos (\omega t)\\ +({\alpha }_1+k_o)c_3\cos (3\omega t)-2n{c}_3\lambda A^{2n-1}{\displaystyle \sum _{s=1}^{\infty }{d}_{2s+1}}\cos ((2s+1)\omega t)\end{array}$$

(36)

Dropping the secular terms in v₂ leads to

$$-{\alpha }_1{c}_3+{\alpha }_{2}A-2n\lambda c_3{d}_1{A}^{2n-1}-k_o{c}_3-k_{o}A+k_{1}A=0$$

(37)

Setting p = 1 in equation (11) and substituting for ${\alpha }_2$ from equation (37) delivers the second-order approximate angular frequency

$${\omega }_2=\sqrt{\alpha +\lambda b_1{A}^{2n-1}+\lambda c_3{A}^{2n-2}(b_1+2n{d}_1)-k_1}$$

(38)

Note that ω₂ does not depend on k_o.

With equation (37) satisfied, equation (36) reduces to

$$v_2^{\prime \prime }+{\omega }^2{v}_2=({\alpha }_1+k_o)c_3\cos (3\omega t)-2n{c}_3\lambda A^{2n-1}{\displaystyle \sum _{s=1}^{\infty }{d}_{2s+1}}\cos ((2s+1)\omega t)$$

(39)

with initial condition $v_2(0)=0$ and $v_2^{\prime }(0)=0$ . The solution of equation (39) can be expressed in the form

$$v_2(t)={\displaystyle \sum _{s=0}^{\infty }{e}_{2s+1}}\cos ((2s+1)\omega t)$$

(40)

Substituting on the left-hand side of equation (39), we get

$$\begin{array}{l}-{\omega }^2{\displaystyle \sum _{s=0}^{\infty }4}s(s+1)e_{2s+1}\cos \left((2s+1)\omega t\right)=({\alpha }_1+k_o)c_3\cos (3\omega t)-2n{c}_3\lambda A^{2n-1}{\displaystyle \sum _{s=1}^{\infty }{d}_{2s+1}}\cos \left(2s+1)\omega t\right)\end{array}$$

(41)

Equating coefficients of corresponding harmonics on both sides we obtain

$$e_3=\frac{(2n{d}_3\lambda A^{2n-1}-{\alpha }_1-k_o)c_3}{8{\omega }^2}$$

(42)

$$e_{2s+1}=\frac{2n{c}_3\lambda A^{2n-1}{d}_{2s+1}}{4s(s+1){\omega }^2},s=1,2,3,\mathrm{\dots }$$

(43)

The initial condition $v_2(0)=0$ delivers

$$e_1=-{\displaystyle \sum _{s=1}^{\infty }{e}_{2s+1}}$$

(44)

Using a similar argument as above we can truncate the series expansion in equation (40)

$$v_2^{(M)}(t)=e_1^{(M)}\cos (\omega t)+{\displaystyle \sum _{s=1}^M{e}_{2s+1}}\cos ((2s+1)\omega t)$$

with $e_1^{(M)}=-{\displaystyle {\sum }_{s=1}^M{e}_{2s+1}}.$

Consider the case M = 2. Then

$$v_2^{(2)}(t)=e_1^{(2)}\cos (\omega t)+e_3\cos (3\omega t)+e_5\cos (5\omega t)$$

(45)

with $e_1^{(2)}=-(e_3+e_5)$ . Substituting above, equation (45) can be written as

$$v_2^{(2)}(t)=e_3(\cos (3\omega t)-\cos (\omega t))+e_5(\cos (5\omega t)-\cos (\omega t))$$

(46)

where e₃ is given by equation (42) with $d_3=\frac{-2(n-7/2)}{(n+5/2)(n+3/2)}{b}_1$ and

$$e_5=\frac{n{c}_3\lambda A^{2n-1}{d}_5}{12{\omega }^2}$$

follows from equation (43) with

$$d_5=-\frac{4n^2-48n+23}{2(n+7/2)(n+5/2)(n+3/2)}{b}_1$$

Finally, the second-order approximate solution of the nonlinear oscillator (equation (8)) is given by

$$\begin{array}{l}{v}_2(t)=v_{°}(t,{\omega }_2)+v_1^{(1)}(t,{\omega }_2)+v_2^{(2)}(t,{\omega }_2)\\ =(A-c_3-e_3-e_5)\cos ({\omega }_{2}t)+(c_3+e_3)\cos (3{\omega }_{2}t)+e_5\cos (5{\omega }_{2}t)\end{array}$$

(47)

where c₃ is given by equation (30) and ${\omega }_2$ is given by equation (38).

Note that ${\omega }_2$ depends nonlinearly on k₁ only. Substituting for ${\alpha }_1$ from equation (22) in equation (42), e₃ does not depend on k_o. Then, the second-order approximate solution (equation (38)) also does not depend on k_o. This leads to an efficient line search method to determine the optimal value of k₁ as already outlined above.

Results and discussion

Consider the nonlinear discontinuity equation (8) with $n=3,\alpha =1$ , and $\lambda =1$ . We take A = 10. Figure 1 shows a plot of the normalized solution ${\overline{v}}_n(t)=v_n(t)/A$ of equation (8) for one period T obtained numerically using the MATLAB ode45 solver. The angular frequency of oscillations was found as ${\omega }_n=\frac{2\pi }{T}=225.2038$ rad s $-1$ .

Figure 1.

Normalized solution of the nonlinear discontinuity equation (8) with $n=3,\alpha =1,\lambda =1$ , and A = 10 for one period obtained numerically using the MATLAB ode45 solver.

The integrated squared error (ISE) of the normalized first-order approximate solution ${\overline{v}}_1(t)=v_1(t)/A$ given by equation (31) was computed as

$${\displaystyle {\int }_0^T({\overline{v}}_1(}t)-{\overline{v}}_n(t){)}^2\mathrm{d}t$$

The ISE with respect to k_o is convex. Figure 2(a) shows a plot of the square root of the ISE for $8000\le k_o\le 8200$ which is minimum at k_o = 8094.

Figure 2.

(a) Square root of the integrated squared error (ISE) for the first-order approximate solution for $8000\le k_o\le 8200$ . The minimizer of the error is k_o = 8094. (b) Absolute value of the error for the first-order approximate solution via standard HPM (k_o = 0) and modified HPM (k_o = 8094).

Substituting k_o = 8094 in equation (23), the first-order approximate angular frequency is $\omega 1=223.8576$ rad s $-1$ with a relative error of 0.6%. In comparison, the standard HPM (k_o = 0) overestimates the angular frequency at 241.2597 rad s $-1$ with a relative error of 7.13%.

Figure 2(b) shows the absolute value of the error over one period for the first-order approximate solution obtained using the modified HPM with k_o = 8094 (dashed line) and the standard HPM with k_o = 0 (solid line). In the case of the modified HPM, the error is bounded below 0.02 whereas for the standard HPM the error is an order of magnitude higher, bounded below 0.33. This result is significant for two reasons: first, the modified HPM with the auxiliary term considerably improves the first-order estimate compared to the standard HPM; second, it emphasizes the need to carefully select the auxiliary parameter. The choice k = 0 ⁴⁷ or k = 1 ⁴⁶ is clearly non-optimal.

Figure 3(a) shows a plot of the square root of the ISE for the second-order approximate solution given by equation (31) for a range of values of the parameter $100\le k_1\le 500$ . We find that the minimizer of the error is $k_1=299$ .

Figure 3.

(a) Square root of the integrated square error (ISE) the second-order approximate solution for $100\le k_1\le 500$ . The minimizer of the error is $k_1=299$ . (b) Absolute value of the error for the second-order approximate solution via standard HPM and modified HPM with $k_1=299$ .

Substituting $k_1=299$ in equation (38), the second-order approximate angular frequency is ${\omega }_2=223.7790$ rad s $-1$ with a relative error of 0.62%. This seems unsatisfactory for two reasons: first, the second-order approximation did not improve on the first-order estimate of the angular frequency; second, the estimate of the angular frequency from the second-order approximation of the standard HPM is 224.5685 rad s $-1$ with a relative error of 0.28%.

However, as shown in Figure 3(b), the absolute value of the error for the second-order approximate solution given by equation (31) with $k_1=299$ is generally better than that obtained using the standard HPM. We find that the error using the modified HPM remains below 0.015 while the error using the standard HPM is considerably higher particularly on the interval $(0.65,1)$ .

Conclusion

In this paper, we have discussed for the first time a systematic procedure for determining the tuning parameter in the homotopy equation of a modified HPM. The optimal value of the parameter is determined as the minimizer of the ISE function. The first- and second-order approximate solutions depend on a single component of the parameter expansion. This leads to an efficient line search method to find the minimizer. The procedure was illustrated for a nonlinear oscillator with a discontinuous term and it was shown that an optimal choice of the auxiliary parameter leads to a significant improvement in the first-order approximation of the angular frequency of oscillations and the solution of the nonlinear equation compared to the standard HPM. While the second-order approximate solution of the angular frequency did not yield an improvement, the second-order approximate solution of the nonlinear equation is generally better. Future work will investigate more general forms of the auxiliary term.