Simplified variational iteration method for solving ordinary differential equations and eigenvalue problems

Abstract

A simplified variational iteration method is proposed to solve high-order homogeneous or nonhomogeneous linear ordinary differential equation and ordinary differential equation eigenvalue problems more efficiently and conveniently. The simplification includes two aspects: (1) explicitly deducing the general form of the differential equation for the identification of the general Lagrange multiplier while avoiding the complexity of variational calculations during identification and (2) simplifying the iterative expressions to reduce the computational work of each iteration. Three ordinary differential equations in mechanics are solved by this simplified variational iteration method, which proves that it is valid and more concise than traditional methods. To make the method more practical, it is suggested that some complicated analytical derivations be executed numerically, thereby achieving a simplified semi-analytical variational iteration method that can be easily implemented by computer programs. The method is then used to numerically solve two complex ordinary differential equation problems derived from the continuum analysis of tall building structures: a sixth-order nonhomogeneous ordinary differential equation with complex boundary conditions and a sixth-order ordinary differential equation eigenvalue problem. Numerical computer programs are developed for these two problems, and corresponding examples are provided to verify the accuracy and efficiency of the simplified variational iteration method in solving complex ordinary differential equation problems.

Introduction

Differential equations are widely used to describe various mechanical problems, ¹ thus making the method used to solve them an important issue in many cases. For low-order and simple differential equations, it is easy to obtain analytical solutions; however, for high-order or complicated ones, analytical solutions are difficult to obtain or may not even exist. Numerical approaches are therefore normally employed in practice to obtain usable results, especially in engineering. ² Ordinary differential equations (ODEs) ³ are equations with only one independent variable and can be mainly divided into two types according to the boundary conditions: initial value problems (IVPs) and boundary value problems (BVPs). Between these two types of problems, the solution of BVPs is more complicated because the boundary conditions are defined at different points. The commonly used numerical solution methods for BVPs are the shooting method, Galerkin method, finite difference method, finite element method, ⁴ and so on. These methods are all capable of solving ODEs, but do have some drawbacks. The shooting method involves solving the BVP’s approximate IVP several times, which could be time-consuming for complicated problems. The Galerkin’s method requires finding a series of trial functions that are compatible with the boundary conditions, but these functions are not always easy to determine. As for the finite difference method and finite element method, the accuracy of the solution is heavily dependent on the density or quality of the mesh. What’s more, these methods are all approximate methods, even for linear ODEs.

Variational iteration is a relatively new method for solving differential equations that is theoretically accurate with simple concepts and fast convergence. It is derived from the general Lagrange multiplier method used for solving nonlinear equations in quantum mechanics, ⁵ which was modified by He^5–8 and Biazar et al. ⁹ into an iteration method named the variational iteration method (VIM). The VIM is in fact a general form of the Newton–Raphson iteration method. ⁵ Moreover, H Jafari ¹⁰ has recently proved that when the highest derivative term is considered the linear part in the VIM while solving nonlinear differential equations, the VIM is equivalent to the well-known classical successive approximations method. The VIM provides an effective means for solving differential equations, large linear systems, and so on and has been used to solve various differential equations in different disciplines, such as the Fokker–Planck equation^11–13 and Boltzmann equation ¹⁴ in statistical mechanics, fractional partial differential equations in fluid mechanics, ¹⁵ nonlinear wave equations,^16–18 partial differential equations for water waves, ¹⁹ variational problems, ²⁰ nonlinear equations in heat transfer,^21,22 free vibration problems of Euler–Bernoulli beams ²³ and marine risers, ²⁴ and eigenvalue problems of Sturm–Liouville differential equations ²⁵ and nonlinear oscillators.^26–28

In this study, the VIM is simplified to solve high-order linear ODEs and ODE eigenvalue problems more conveniently and practically. First, the general form of the differential equation for identifying the general Lagrange multiplier is explicitly given, and then the iterative formulae are simplified for linear ODEs. In addition, integral expressions are suggested to be performed numerically in order to avoid complicated analytical derivations and make the method more convenient for numerical computer programming. Finally, several ODEs are selected for solving by the simplified VIM to verify its validity and efficiency.

Simplification of VIM

In general, an ODE can be expressed in the following form

$$\begin{array}{cc}F[y\left(z\right),z]=0\hfill & z\in [a,b]\hfill \end{array}$$

(1)

where F denotes the differential operator of the ODE.

Now, suppose that the solution of the above ODE can be expressed in an iterative form as

$$y_n\left(z\right)=y_{n-1}\left(z\right)+\underset{a}{\overset{z}{\int }}\sigma (x,z)·F[y_{n-1}\left(x\right),x]\mathrm{d}x$$

(2)

where $\sigma (x,z)$ is the generalized Lagrange multiplier and $y_n(z)\mathrm{and}{y}_{n-1}(z)$ denote the solution after n and n−1 iterations, respectively.

The real solution is denoted by $y(z)$ (so that $F[y(z)]=0$ ), and the discrepancy between $y_{n-1}(z)$ and $y(z)$ is $\delta y_{n-1}(z)$ , that is

$$y_{n-1}=y+\delta y_{n-1}$$

(3)

Substituting the above equation in equation (2), we have

$$y_n=y+\delta y_{n-1}+\underset{a}{\overset{z}{\int }}\sigma ·F(y+\delta y_{n-1})\mathrm{d}x$$

(4)

Supposing that $\delta y_{n-1}$ is sufficiently small, the first-order Taylor series expansion of $F(y+\delta y_{n-1})$ gives

$$\begin{array}{cc}\multicolumn{1}{c}{y_n}& =y+\delta y_{n-1}+\underset{a}{\overset{z}{\int }}\sigma [F\left(y\right)+\delta y_{n-1}\frac{\partial F\left(y\right)}{\partial y}]\mathrm{d}\epsilon +O\left(\delta y_{n-1}^2\right)\hfill \\ \multicolumn{1}{c}{}& =y+\delta y_{n-1}+\underset{a}{\overset{z}{\int }}\sigma \left[\delta y_{n-1}\frac{\partial F\left(y\right)}{\partial y}\right]\mathrm{d}\epsilon +O\left(\delta y_{n-1}^2\right)\hfill \end{array}$$

(5)

Provided that

$$\delta y_{n-1}+\underset{a}{\overset{z}{\int }}\sigma \left[\delta y_{n-1}\frac{\partial F\left(y\right)}{\partial y}\right]\mathrm{d}x=0$$

(6)

holds, the discrepancy after n iterations $\delta y_n$ is a second-order small quantity of the discrepancy after n−1 iterations $\delta y_{n-1}$ . This means that the nth iterative solution $y_n$ is closer to the real solution than the (n−1)th iterative solution $y_{n-1}$ .

Equation (6) can therefore be used to identify the optimum generalized Lagrange multiplier $\sigma (x,z)$ for the iterative solution. However, as the equation includes the unknown real solution $y(z)$ that makes it unusable, an approximation must be made herein for practice. That is, $y(z)$ is replaced with $y_{n-1}(z)$

$$\delta y_{n-1}+\underset{a}{\overset{z}{\int }}\sigma ·\delta F\left(y_{n-1}\right)\mathrm{d}x=0$$

(7)

If we now consider a special case in which $F(y)$ is not a functional but a function, then equation (7) can be simplified as

$$\mathrm{d}{y}_{n-1}+\Lambda F\prime \left(y_{n-1}\right)\mathrm{d}{y}_{n-1}=0\Rightarrow \Lambda =-{[F\prime \left(y_{n-1}\right)]}^{-1}$$

(8)

That is

$$y_n=y_{n-1}+\Lambda F\left(y_{n-1}\right)=y_{n-1}-\frac{F\left(y_{n-1}\right)}{F\prime \left(y_{n-1}\right)}$$

(9)

The above formula is the iteration formula of the famous Newton–Raphson method for solving nonlinear equations, which means that the VIM can be regarded as a general form of the Newton–Raphson method or that the Newton–Raphson method is a particular case of the VIM. In this manner, the validity and feasibility of the VIM is also verified from a unique perspective.

When the operator F is complicated, it is difficult to obtain the expression for $\sigma (x,z)$ using equation (6), and so some approximations need to be made. To begin with, the operator F is divided into three parts

$$F[y\left(z\right),z]=L\left[y\left(z\right)\right]+N\left[\tilde{y}\left(z\right)\right]+G\left(z\right)$$

(10)

where L is a simple linear operator (so simple that the analytical solution of $L[y(z)]=0$ is easy to obtain); N is a complicated linear or nonlinear operator; $\tilde{y}$ is the restricted variation, that is, $\delta \tilde{y}=0$ , such that $\delta N[\tilde{y}(z)]=0$ ; and $G(z)$ is a function of z and independent of y, that is, $\delta G=0$ . Therefore, equation (7) can be rewritten as

$$\delta y_{n-1}\left(z\right)+\underset{a}{\overset{z}{\int }}\sigma (x,z)·\delta L\left[y_{n-1}\left(x\right)\right]\mathrm{d}x=0$$

(11)

Equation (11) produces several stationary conditions that can be treated as a differential equation (equation (12)) together with the corresponding boundary conditions (equation (13)) with respect to the generalized Lagrange multiplier $\sigma (x,z)$

$$L\left[\sigma \left(x\right)\right]=0$$

(12)

$$\begin{array}{c}1-L_{-1}{\left[\sigma \left(x\right)\right]}_{x=z}=0\hfill \\ L_{-2}{\left[\sigma \left(x\right)\right]}_{x=z}=0\hfill \\ L_{-3}{\left[\sigma \left(x\right)\right]}_{x=z}=0\hfill \\ \cdots \cdots \hfill \\ L_{-(n-1)}{\left[\sigma \left(x\right)\right]}_{x=z}=0\hfill \\ {\sigma \left(x\right)|}_{x=z}=0\hfill \end{array}$$

(13)

where $L_{-k}[•]$ indicates that the differential order of all terms is reduced by k in the operator L and n is the highest order in the operator L; if $k>n$ , then let $L_{-k}[•]=0$ . In fact, equations (12) and (13) are equivalent to the method suggested by Jafari et al. ²⁹

Later, $\sigma (x,z)$ will be identified by solving the differential equation expressed in equation (12), and then, equation (2) can be used to perform the iterations for the solution of equation (1). The solution of $L[y(z)]=0$ can be chosen as the initial solution $y_0$ ; if the iteration count is larger than one, the boundary conditions at $z=a$ (the starting point of the integral) should be considered when solving $L[y(z)]=0$ because the iterations will not change the boundary value of y at $z=a$ in equation (2). Neglecting the operator N means that the resulting generalized Lagrange multiplier $\sigma (x,z)$ is not the optimal one, but rather a practical one.

Based on the above derivations, the operator L and Lagrange multiplier $\sigma (x,z)$ possess the following properties

$$\begin{array}{cc}L\left[y_0\left(z\right)\right]=0\hfill & L\left[\sigma (x,z)\right]=0\hfill \end{array}$$

(14)

$${\frac{{\partial }^m\sigma (x,z)}{\partial x^m}|}_{x=z}=\{\begin{array}{cc}0\hfill & 0\le m<n-1\hfill \\ \multicolumn{1}{c}{-1}& m=n-1\hfill \end{array}$$

(15)

In view of the derivation formula for the uncertain limit integral given below

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\mathrm{d}}{\mathrm{d}z}\left[\underset{\phi \left(z\right)}{\overset{\psi \left(z\right)}{\int }}f(x,z)\mathrm{d}x\right]=f[z,\psi \left(z\right)]\frac{\mathrm{d}\psi \left(z\right)}{\mathrm{d}z}}\\ \multicolumn{1}{c}{-f[z,\phi \left(z\right)]\frac{\mathrm{d}\phi \left(z\right)}{\mathrm{d}z}+\underset{\phi \left(z\right)}{\overset{\psi \left(z\right)}{\int }}\frac{\partial f(x,z)}{\partial z}\mathrm{d}x}\end{array}$$

(16)

We can obtain a new property of $\sigma (x,z)$

$$\frac{{\mathrm{d}}^m}{\mathrm{d}{z}^m}{\displaystyle \underset{a}{\overset{z}{\int }}\sigma \left(x,z\right)\cdot r\left(x\right)\mathrm{d}x}=\{\begin{array}{cc}{\displaystyle \underset{a}{\overset{z}{\int }}\frac{{\partial }^m\sigma \left(x,z\right)}{\partial z^m}\cdot r\left(x\right)\mathrm{d}x}& 0<m<n\\ -r\left(z\right)+{\displaystyle \underset{a}{\overset{z}{\int }}\frac{{\partial }^m\sigma \left(x,z\right)}{\partial z^m}\cdot r\left(x\right)\mathrm{d}x}& m=n\end{array}$$

(17)

Because the linear operator L can be expanded as

$$L\left[y\left(z\right)\right]=\sum _{m=0}^n{c}_m\frac{{\mathrm{d}}^{m}y}{\mathrm{d}{z}^m}$$

(18)

We obtain

$$\begin{array}{c}L[\underset{a}{\overset{z}{\int }}\sigma (x,z)·r\left(x\right)\mathrm{d}x]=\sum _{m=0}^n{c}_m\frac{{\mathrm{d}}^m}{\mathrm{d}{z}^m}\underset{a}{\overset{z}{\int }}\sigma (x,z)·r\left(x\right)\mathrm{d}x\hfill \\ =-r\left(z\right)+\underset{a}{\overset{z}{\int }}\sum _{m=0}^n{c}_m\frac{{\partial }^m\sigma (x,z)}{\partial z^m}·r\left(x\right)\mathrm{d}x\hfill \\ =-r\left(z\right)+\underset{a}{\overset{z}{\int }}L\left[\sigma (x,z)\right]·r\left(x\right)\mathrm{d}x\hfill \\ =-r\left(z\right)\hfill \end{array}$$

(19)

That is

$$L[\underset{a}{\overset{z}{\int }}\sigma (x,z)·r\left(x\right)\mathrm{d}x]+r\left(z\right)=0$$

(20)

We now define a new symbol ${\nu }_n(z)$ as

$${\nu }_n\left(z\right)=\underset{a}{\overset{z}{\int }}\sigma (x,z)·\left\{N\left[{\nu }_{n-1}\left(x\right)\right]\right\}\mathrm{d}x$$

(21)

where ${\nu }_0(z)=y_0(z)$ . Substituting $r(z)={\nu }_n(z)$ in equation (20), the following equation can be obtained

$$L\left[{\nu }_n\left(z\right)\right]+N\left[{\nu }_{n-1}\left(z\right)\right]=0$$

(22)

Supposing that N is also a linear operator, as per equation (2), we get

$$\begin{array}{c}{y}_1\left(z\right)=y_0\left(z\right)+{\displaystyle \underset{a}{\overset{z}{\int }}\sigma \left(x,z\right)\cdot \left\{L\left[y_0\left(x\right)\right]+N\left[y_0\left(x\right)\right]+G\left(x\right)\right\}\mathrm{d}x}\\ =y_0\left(z\right)+{\displaystyle \underset{a}{\overset{z}{\int }}\sigma \left(x,z\right)\cdot N\left[y_0\left(x\right)\right]\mathrm{d}x}+{\displaystyle \underset{a}{\overset{z}{\int }}\sigma \left(x,z\right)\cdot G\left(x\right)\mathrm{d}x}\\ ={\nu }_0\left(z\right)+{\nu }_1\left(z\right)+{\displaystyle \underset{a}{\overset{z}{\int }}\sigma \left(x,z\right)\cdot G\left(x\right)\mathrm{d}x}\\ y_2\left(z\right)=y_1\left(z\right)+{\displaystyle \underset{a}{\overset{z}{\int }}\sigma \left(x,z\right)\cdot \left\{L\left[y_1\left(x\right)\right]+N\left[y_1\left(x\right)\right]+G\left(x\right)\right\}\mathrm{d}x}\\ =y_1\left(z\right)+{\displaystyle \underset{a}{\overset{z}{\int }}\sigma \left(x,z\right)\cdot \left\{\begin{array}{l}L\left[{\nu }_0\left(x\right)\right]+L\left[{\nu }_1\left(x\right)\right]\hfill \\ +N\left[{\nu }_0\left(x\right)\right]+N\left[{\nu }_1\left(x\right)\right]+G\left(x\right)\hfill \end{array}\right\}\mathrm{d}x}\\ =y_1\left(z\right)+{\displaystyle \underset{a}{\overset{z}{\int }}\sigma \left(x,z\right)\cdot \left\{N\left[{\nu }_1\left(x\right)\right]+G\left(x\right)\right\}\mathrm{d}x}\\ ={\nu }_0\left(z\right)+{\nu }_1\left(z\right)+{\nu }_2\left(z\right)+{\displaystyle \underset{a}{\overset{z}{\int }}\sigma \left(x,z\right)\cdot G\left(x\right)\mathrm{d}x}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{y_3\left(z\right)=y_2\left(z\right)+\underset{a}{\overset{z}{\int }}\sigma (x,z)·\{L\left[y_2\left(x\right)\right]+N\left[y_2\left(x\right)\right]+G\left(x\right)\}\mathrm{d}x}\\ =y_1\left(z\right)+\underset{a}{\overset{z}{\int }}\sigma (x,z)·\left\{\begin{array}{c}L\left[{\nu }_0\left(x\right)\right]+L\left[{\nu }_1\left(x\right)\right]+L\left[{\nu }_2\left(x\right)\right]\hfill \\ \multicolumn{1}{c}{\begin{array}{c}+N\left[{\nu }_0\left(x\right)\right]+N\left[{\nu }_1\left(x\right)\right]\hfill \\ \multicolumn{1}{c}{+N\left[{\nu }_2\left(x\right)\right]+G\left(x\right)}\end{array}}\end{array}\right\}\mathrm{d}x\hfill \\ =y_1\left(z\right)+\underset{a}{\overset{z}{\int }}\sigma (x,z)·\{N\left[{\nu }_2\left(x\right)\right]+G\left(x\right)\}\mathrm{d}x\hfill \\ ={\nu }_0\left(z\right)+{\nu }_1\left(z\right)+{\nu }_2\left(z\right)+{\nu }_3\left(z\right)+\underset{a}{\overset{z}{\int }}\sigma (x,z)·G\left(x\right)\mathrm{d}x\hfill \\ \multicolumn{1}{c}{\cdots }\end{array}$$

Analogically, the ith iterative solution can be expressed as

$$y_i\left(z\right)=\sum _{k=0}^i{\nu }_k\left(z\right)+\underset{a}{\overset{z}{\int }}\sigma (x,z)·G\left(x\right)\mathrm{d}x$$

(24)

The above expression means that iterations need not be carried out using equation (2), but instead only ${\nu }_1(z),{\nu }_2(z),{\nu }_3(z)\cdots {\nu }_i(z)$ needs to be calculated in accordance with equation (21). This will eliminate the operator L, making the iterative solving procedure simpler and easier.

From equations (21) and (24), we can also infer that when $N[y(z)]=0$ , ${\nu }_k(z)$ also equals 0. An accurate analytical solution can then be obtained after only one iteration

$$y\left(z\right)=y_0\left(z\right)+\underset{a}{\overset{z}{\int }}\sigma (x,z)·G\left(x\right)\mathrm{d}x$$

(25)

For IVPs, the initial solution $y_0(z)$ and every solution obtained during the iterative solving procedure can be fully confirmed according to the initial conditions without undetermined coefficients. Thus, if the iteration converges, the solution calculated by equation (24) can be used directly. With BVPs, however, the initial solution $y_0(z)$ contains undetermined coefficients. When the iteration is finished, these undetermined coefficients need to be eliminated using the remaining boundary conditions. The initial solution with undetermined coefficients can generally be written as

$$y_0\left(z\right)=\sum _{k=1}^q{C}_k{f}_{k,0}\left(z\right)$$

(26)

where $C_k$ denotes the undetermined coefficient, $f_{k,0}(z)$ denotes one of the particular solutions of $L[y(z)]=0$ , and q is the number of particular solutions.

The ith iterative solution $y_i(z)$ can be rewritten as follows (the derivation is similar to that of equation (23))

$$y_i\left(z\right)=\sum _{k=1}^q{C}_k\left[\sum _{l=0}^i{f}_{k,l}\left(z\right)\right]+\underset{a}{\overset{z}{\int }}\sigma (x,z)·G\left(x\right)\mathrm{d}x$$

(27)

where

$$f_{k,l}\left(z\right)=\underset{a}{\overset{z}{\int }}\sigma (x,z)·N\left[f_{k,l-1}\left(x\right)\right]\mathrm{d}x$$

(28)

Equation (27) divides the iterative solution into several components, each of which can be calculated independently (equation (28)), so the undetermined coefficients can be ignored during the iterations. In this manner, the derivation in every iterative step is kept brief and simple, making it easy to obtain the numerical approaches involved.

As for the eigenvalue problems of ODEs, the operator N can usually be written in the following form

$$N\left[y\left(z\right)\right]=\lambda \overline{N}\left[y\left(z\right)\right]$$

(29)

By repeating the deduction in equation (23), the solution after the ith iteration can be obtained as

$$y_i\left(z\right)=\sum _{k=1}^q{C}_k\left[\sum _{l=0}^i{\lambda }^l{f}_{k,l}\left(z\right)\right]+\underset{a}{\overset{z}{\int }}\sigma (x,z)·G\left(x\right)\mathrm{d}x$$

(30)

where

$$f_{k,l}\left(z\right)=\underset{a}{\overset{z}{\int }}\sigma (x,z)·\overline{N}\left[f_{k,l-1}\left(x\right)\right]\mathrm{d}x$$

(31)

Verification of simplified VIM

Here, three simple but classical problems will be solved using the simplified VIM to demonstrate the solution procedure and verify the validity of the simplification.

Problem I: Forced vibration of undamped single-degree-of-freedom system

The differential equation of forced vibration for an undamped single-degree-of-freedom (SDOF) system is ²⁷

$$\stackrel{··}{u}\left(t\right)+{\omega }^{2}u\left(t\right)=\frac{1}{m}p\left(t\right)$$

(32)

whose initial conditions are

$$\begin{array}{cc}u\left(0\right)=d_0\hfill & \stackrel{·}{u}\left(0\right)=v_0\hfill \end{array}$$

(33)

According to equation (10), the operators can be identified as

$$\begin{array}{ccc}L\left[u\left(t\right)\right]=\stackrel{··}{u}\left(t\right)+{\omega }^{2}u\left(t\right)\hfill & N\left[u\left(t\right)\right]=0\hfill & G\left(t\right)=-\frac{1}{m}p\left(t\right)\hfill \end{array}$$

(34)

Substituting the above equations in equations (12) and (13), we have

$$\begin{array}{c}\sigma ″\left(x\right)+{\omega }^2\sigma \left(x\right)=0\hfill \\ \begin{array}{cc}1-\sigma \prime \left(t\right)=0\hfill & \sigma \left(t\right)=0\hfill \end{array}\hfill \end{array}$$

(35)

The generalized Lagrange multiplier $\sigma (x,z)$ can then be solved

$$\sigma (x,t)=\frac{\sin \left[\omega (x-t)\right]}{\omega }$$

(36)

Forcing $L[u(t)]=0$ to obtain the initial solution (considering the conditions in equation (33)), we get

$$u_0\left(t\right)=d_0\cos \omega t+\frac{v_0\sin \omega t}{\omega }$$

(37)

Because $N[u(t)]=0$ , the accurate analytical solution can be determined according to equation (25)

$$\begin{array}{l}u\left(t\right)=d_0\cos \omega t+\frac{v_0\sin \omega t}{\omega }\\ +\frac{1}{m\omega }{\displaystyle \underset{0}{\overset{t}{\int }}\sin \left[\omega \left(t-x\right)\right]·p\left(x\right)\mathrm{d}x}\end{array}$$

(38)

Equation (38) is the well-known Duhamel integral or the solution of the vibration of an SDOF system under arbitrary dynamic loads. This proves that the identified generalized Lagrange multiplier based on equations (12) and (13) is rational, and that the simplified VIM is valid and capable of solving IVPs.

Problem II: Cantilever (Euler–Bernoulli beam) under distributed transverse load

If the distributed transverse load is set to $p(\zeta )=\mathrm{EI}/H^4\sin (\pi \zeta /2)$ , then the static equilibrium differential equation of a cantilever is

$$\frac{{\mathrm{d}}^{4}y}{\mathrm{d}{\zeta }^4}-\sin \left(\frac{\pi \zeta }{2}\right)=0$$

(39)

where $\zeta \in [0,1]$ . The corresponding boundary conditions are

$$\begin{array}{cccc}y\left(0\right)=0\hfill & y\prime \left(0\right)=0\hfill & y″\left(1\right)=0\hfill & y″\prime \left(1\right)=0\hfill \end{array}$$

(40)

The analytical solution of equation (39) can be obtained by simply integrating the equation four times

$$y\left(\zeta \right)=\frac{2}{{\pi }^4}[{\pi }^2{\zeta }^2-4\pi \zeta +8\sin \left(\frac{\pi \zeta }{2}\right)]$$

(41)

To now try and solve equation (39) using the simplified VIM, the operators should first be identified

$$L\left(y\right)=y^{\left(4\right)},N\left(y\right)=0,G\left(\zeta \right)=-\sin (\pi \zeta /2)$$

(42)

Then, based on equations (12) and (13), we have

$$\begin{array}{c}{\sigma }^{(4)}\left(x\right)=0\hfill \\ \begin{array}{cccc}1-\sigma ″\prime \left(\zeta \right)=0\hfill & \sigma ″\left(\zeta \right)=0\hfill & \sigma \prime \left(\zeta \right)=0\hfill & \sigma \left(\zeta \right)=0\hfill \end{array}\hfill \end{array}$$

(43)

which yields the generalized Lagrange multiplier

$$\sigma (x,\zeta )=\frac{{(x-\zeta )}^3}{6}$$

(44)

Considering $L[y(\zeta )]=0$ , the initial solution can be found

$$y_0=C_1+C_2\zeta +C_3{\zeta }^2+C_4{\zeta }^3$$

(45)

Because $N[y(\zeta )]=0$ , the final solution can be determined according to equation (25) as

$$\begin{array}{cc}\multicolumn{1}{c}{y\left(\zeta \right)}& =y_0\left(\zeta \right)-\underset{0}{\overset{\zeta }{\int }}\frac{{(x-\zeta )}^3}{6}·\sin \left(\frac{\pi \zeta }{2}\right)\mathrm{d}x\hfill \\ \multicolumn{1}{c}{}& =C_1+C_2\zeta +C_3{\zeta }^2+C_4{\zeta }^3\hfill \\ \multicolumn{1}{c}{}& +\frac{1}{3{\pi }^4}[{\pi }^3{\zeta }^3-24\pi \zeta +48\sin (\pi \zeta /2)]\hfill \end{array}$$

(46)

The undetermined coefficients can be obtained using the boundary conditions in equation (40), as per the following

$$C_1=C_2=0,C_3=\frac{2}{{\pi }^2},C_4=-\frac{1}{3\pi }$$

(47)

If we substitute the value of the above coefficients in equation (46), the expression obtained will be the same as that in equation (41). This example proves that the simplified VIM is capable of solving BVPs.

Problem III: Free vibration of cantilever (Euler–Bernoulli beam)

The eigenvalue differential equation of transverse free vibration of a cantilever (Euler–Bernoulli beam) is

$$y^{\left(4\right)}\left(\zeta \right)-{\left(\lambda H\right)}^{4}y\left(\zeta \right)=0$$

(48)

The boundary conditions are the same as those in equation (40). The precise results of eigenvalues are the roots of the following equation

$$1+\cos \left(\lambda H\right)\cosh \left(\lambda H\right)=0$$

(49)

The numerical values of the accurate eigenvalues (the first two orders) are

$${\left(\lambda H\right)}_1=1.8751,{\left(\lambda H\right)}_2=4.6941$$

(50)

When using the VIM, the operator L is chosen as $L(y)=y^{(4)}$ . As this is the same as that used in the previous problem, the generalized Lagrange multiplier will also be the same as the one in equation (44). In addition, the initial solution should be the same as the one in equation (45). To reduce the computational burden, the initial solution is simplified in view of the boundary conditions at $\zeta =0$ in equation (40)

$$y_0=C_3{\zeta }^2+C_4{\zeta }^3$$

(51)

The other operators are can then be identified as

$$N\left[y\left(\zeta \right)\right]=-{\left(\lambda H\right)}^{4}y\left(\zeta \right),G\left(\zeta \right)=0$$

(52)

Because $L\mathrm{and}N$ are both linear operators, equation (28) can be applied here

$$\begin{array}{cc}{f}_{1,0}={\zeta }^2\hfill & f_{2,0}={\zeta }^3\hfill \\ \multicolumn{1}{c}{f_{1,1}=\frac{{\left(\lambda H\right)}^4}{360}{\zeta }^6}& f_{2,1}=\frac{{\left(\lambda H\right)}^4}{840}{\zeta }^7\hfill \\ \multicolumn{1}{c}{f_{2,1}=\frac{{\left(\lambda H\right)}^8{\zeta }^{10}}{1814400}}& f_{2,2}=\frac{{\left(\lambda H\right)}^8{\zeta }^{11}}{6652800}\hfill \\ \multicolumn{1}{c}{f_{3,1}=\frac{{\left(\lambda H\right)}^{12}{\zeta }^{14}}{43589145600}}& f_{3,2}=\frac{{\left(\lambda H\right)}^{12}{\zeta }^{15}}{217945728000}\hfill \\ \multicolumn{1}{c}{\cdots }& \cdots \hfill \end{array}$$

(53)

Substituting the above equations, we get

$$\begin{array}{c}{y}_1\left(\zeta \right)=C_3·[{\zeta }^2+\frac{{\left(\lambda H\right)}^4}{360}{\zeta }^6]+C_4·[{\zeta }^3+\frac{{\left(\lambda H\right)}^4}{840}{\zeta }^7]\hfill \\ y_2\left(\zeta \right)=C_3·[{\zeta }^2+\frac{{\left(\lambda H\right)}^4{\zeta }^6}{360}+\frac{{\left(\lambda H\right)}^8{\zeta }^{10}}{1814400}]\hfill \\ +C_4·[{\zeta }^3+\frac{{\left(\lambda H\right)}^4{\zeta }^7}{840}+\frac{{\left(\lambda H\right)}^8{\zeta }^{11}}{6652800}]\hfill \\ y_3\left(\zeta \right)=C_3·[{\zeta }^2+\frac{{\left(\lambda H\right)}^4{\zeta }^6}{360}+\frac{{\left(\lambda H\right)}^8{\zeta }^{10}}{1814400}+\frac{{\left(\lambda H\right)}^{12}{\zeta }^{14}}{43589145600}]\hfill \\ +C_4·[{\zeta }^3+\frac{{\left(\lambda H\right)}^4{\zeta }^7}{840}+\frac{{\left(\lambda H\right)}^8{\zeta }^{11}}{6652800}+\frac{{\left(\lambda H\right)}^{12}{\zeta }^{15}}{217945728000}]\hfill \\ \multicolumn{1}{c}{\cdots }\end{array}$$

(54)

Considering the boundary conditions at $z=1$ in equation (40), a polynomial equation of $\lambda H$ is determined (by considering only the first three iterative solutions)

$$\begin{array}{c}{\left(\lambda H\right)}^8-240{\left(\lambda H\right)}^4+2880=0\hfill \\ {\left(\lambda H\right)}^{16}-672{\left(\lambda H\right)}^{12}+2903040{\left(\lambda H\right)}^8\hfill \\ -1219276800{\left(\lambda H\right)}^4+14631321600=0\hfill \\ {\left(\lambda H\right)}^{24}+7920{\left(\lambda H\right)}^{20}+20243520{\left(\lambda H\right)}^{16}\hfill \\ -199264665600{\left(\lambda H\right)}^{12}\hfill \\ +591816056832000{\left(\lambda H\right)}^8-248562743869440000{\left(\lambda H\right)}^4\hfill \\ +2982752926433280000=0\hfill \end{array}$$

(55)

The numerical roots of the above equations (Table 1) indicate that the iterative result will converge step-by-step to an accurate solution as the iteration count increases. Thus, although the simplified VIM can also solve the eigenvalues of ODEs, the expression of the solution will become increasingly complicated (i.e. the number of terms may grow exponentially) as the iterations proceed. To make the application of the method more convenient, we suggest that the integrals and differentials in equation (28) be calculated by numerical methods instead of analytical derivation, which would make the method a semi-analytical but practical one. In addition, numerical differentials should be avoided because these make it difficult to guarantee accuracy.

Table 1.

Eigenvalues of free vibration of cantilever (Euler–Bernoulli beam) calculated by simplified variational iteration method.

Iteration	First-order eigenvalue (accurate value: 1.8751)		Second-order eigenvalue (accurate value: 4.6941)
	Iterative result	Relative error (%)	Iterative result	Relative error (%)
1	1.8866	0.6	3.8830	17.3
2	1.8751	0.0	4.5363	3.4
3	1.8751	0.0	4.6897	0.1

Application of simplified VIM to complicated ODE problems

In this section, two complicated ODE problems are solved by the proposed VIM method in order to verify and investigate the feasibility and convenience of the method.

Problem description

In the classic theory of tall building analysis, structures are usually idealized as a continuous sandwich cantilevers ³⁰ so that their static and dynamic behavior can be described using differential equations.

The following differential equations are those used to describe a tall building structure under static load and free vibration, respectively ³⁰

$$\frac{{\mathrm{d}}^{6}w}{\mathrm{d}{\zeta }^6}-{\left(k\alpha H\right)}^2\frac{{\mathrm{d}}^{4}w}{\mathrm{d}{\zeta }^4}=\frac{H^4}{\mathrm{EI}}[\frac{{\mathrm{d}}^{2}p}{\mathrm{d}{\zeta }^2}-{\left(k\alpha H\right)}^2(1-k^{-2})p\left(\zeta \right)]$$

(56)

$$\begin{array}{l}\frac{{\partial }^{6}w}{\partial {\zeta }^6}-{\left(k\alpha H\right)}^2\frac{{\partial }^{4}w}{\partial {\zeta }^4}+\frac{m{H}^4}{EI}\\ \left[\frac{{\partial }^{4}w}{\partial {\zeta }^2\partial t^2}-{\left(k\alpha H\right)}^2\left(1-k^{-2}\right)\frac{{\partial }^{2}w}{\partial t^2}\right]=0\end{array}$$

(57)

where, $k,\alpha ,H,m,\mathrm{and}p$ are parameters of a tall building, $\zeta$ is the independent variable, w is the undetermined function, and t is the time-domain variable in dynamic analysis.

Equation (57) can be changed into an ODE by spectral transformation in order to solve the free vibration eigenvalues of a tall building structure

$$\frac{{\mathrm{d}}^{6}y}{\mathrm{d}{\zeta }^6}-{\left(k\alpha H\right)}^2\frac{{\mathrm{d}}^{4}y}{\mathrm{d}{\zeta }^4}-{\left(\lambda H\right)}^4[\frac{{\mathrm{d}}^{2}y}{\mathrm{d}{\zeta }^2}-{\left(k\alpha H\right)}^2(1-k^{-2})y]=0$$

(58)

where $\lambda H$ is the free vibration eigenvalue.

The corresponding boundary conditions of the ODE (56) are

$$\begin{array}{c}\begin{array}{cc}{w|}_{\zeta =0}=0\hfill & \frac{1}{{\left(\beta H\right)}^2}{·\frac{{\mathrm{d}}^{2}w}{\mathrm{d}{\zeta }^2}|}_{\zeta =0}-{\frac{\mathrm{d}w}{\mathrm{d}\zeta }|}_{\zeta =0}=0\hfill \end{array}\hfill \\ {\frac{{\mathrm{d}}^{3}w}{\mathrm{d}{\zeta }^3}|}_{\zeta =0}-\frac{{\left(k\alpha H\right)}^2}{k^2}{\frac{\mathrm{d}w}{\mathrm{d}\zeta }|}_{\zeta =0}+\frac{H^4}{\mathrm{EI}}\underset{0}{\overset{1}{\int }}p\left(\zeta \right)\mathrm{d}\zeta =0\hfill \\ \begin{array}{cc}{\frac{{\mathrm{d}}^{2}w}{\mathrm{d}{\zeta }^2}|}_{\zeta =1}=0\hfill & {\frac{{\mathrm{d}}^{4}w}{\mathrm{d}{\zeta }^4}|}_{\zeta =1}-{\frac{p{H}^4}{\mathrm{EI}}|}_{\zeta =1}=0\hfill \end{array}\hfill \\ {\frac{{\mathrm{d}}^{5}w}{\mathrm{d}{\zeta }^5}|}_{\zeta =1}-{\left(k\alpha H\right)}^2{\frac{{\mathrm{d}}^{3}w}{\mathrm{d}{\zeta }^3}|}_{\zeta =1}\hfill \\ -\frac{H^4}{\mathrm{EI}}[{\frac{\mathrm{d}p}{\mathrm{d}\zeta }|}_{\zeta =1}+{\left(k\alpha H\right)}^2(1-k^{-2})\frac{P}{H}]=0\hfill \end{array}$$

(59)

and these are the boundary conditions for the ODE (58)

$$\begin{array}{c}\begin{array}{cc}{y|}_{\zeta =0}=0\hfill & \frac{1}{{\left(\beta H\right)}^2}·{\frac{{\mathrm{d}}^{2}y}{\mathrm{d}{\zeta }^2}|}_{\zeta =0}-{\frac{\mathrm{d}y}{\mathrm{d}\zeta }|}_{\zeta =0}=0\hfill \end{array}\hfill \\ {\frac{{\mathrm{d}}^{3}y}{\mathrm{d}{\zeta }^3}|}_{\zeta =0}-\frac{{\left(k\alpha H\right)}^2}{k^2}{\frac{\mathrm{d}y}{\mathrm{d}\zeta }|}_{\zeta =0}+{\left(\lambda H\right)}^4\underset{0}{\overset{1}{\int }}y\mathrm{d}\zeta =0\hfill \\ \begin{array}{cc}{\frac{{\mathrm{d}}^{2}y}{\mathrm{d}{\zeta }^2}|}_{\zeta =1}=0\hfill & {\frac{{\mathrm{d}}^{4}y}{\mathrm{d}{\zeta }^4}|}_{\zeta =1}-{\left(\lambda H\right)}^4{y|}_{\zeta =1}=0\hfill \end{array}\hfill \\ {\frac{{\mathrm{d}}^{5}y}{\mathrm{d}{\zeta }^5}|}_{\zeta =1}-{\left(k\alpha H\right)}^2{\frac{{\mathrm{d}}^{3}y}{\mathrm{d}{\zeta }^3}|}_{\zeta =1}-{\left(\lambda H\right)}^4{\frac{\mathrm{d}y}{\mathrm{d}\zeta }|}_{\zeta =1}=0\hfill \end{array}$$

(60)

The general procedure used to obtain the analytical solution of the eigenvalue differential equation (equation (58)) is as follows. First, treat the equation as a normal homogeneous ODE with constant coefficients and find its general solution (include six integration constants $C_1,C_2,\dots ,C_6$ and the eigenvalue $\lambda H$ in the expression). Then, substitute the general solution in equation (60) to obtain a system of linear equations of integration constants $C_1,C_2,\dots ,C_6$ . Next, obtain the equation for $\lambda H$ (so-called eigenvalue equation) by applying Cramer’s rule and extract eigenvalues by solving the equation. Although this procedure seems feasible in theory, it is difficult in practice because the characteristic equation (equation (58)) is a univariate sextic equation (which can be reduced to be quintic), which means that the roots of the equation cannot be expressed in algebraic form according to Abel’s impossibility theorem. In addition, even if a computer algebra system is employed to acquire an expression for the general solution of equation (58), such acquisition will be very cumbersome. It will also include complicated operations on the eigenvalue, $\lambda H$ , such as exponentiation and rooting, which will make it difficult to obtain the value of $\lambda H$ even when using a numerical approach. In view of the difficulty involved in solving equation (58), and the shortcomings of traditional methods mentioned in the Introduction section, the semi-analytical method based on the VIM deduced from previous discussions is used to solve the problem. As for equation (56), although the analytical solution is not very difficult to obtain if the lateral load is analytical, the semi-analytical method will still be applied for convenience.

Solution of sixth-order nonhomogeneous ODE

As mentioned above, the first task involved in solving differential equations by the VIM is to identify the generalized Lagrange multiplier. Before performing this task, the operators in equation (10) should be determined in comparison with equation (56) as follows

$$\begin{array}{c}\begin{array}{cc}L\left[w\left(\zeta \right)\right]=\frac{{\mathrm{d}}^{6}w}{\mathrm{d}{\zeta }^6}-{\left(k\alpha H\right)}^2\frac{{\mathrm{d}}^{4}w}{\mathrm{d}{\zeta }^4}\hfill & N\left[w\left(\zeta \right)\right]=0\hfill \end{array}\hfill \\ G\left(\zeta \right)=\frac{H^4}{\mathrm{EI}}[\frac{{\mathrm{d}}^{2}p\left(\zeta \right)}{\mathrm{d}{\zeta }^2}-{\left(k\alpha H\right)}^2(1-k^{-2})p\left(\zeta \right)]\hfill \end{array}$$

(61)

Then, the recognition equation of the generalized Lagrange multiplier can be obtained based on equations (12) and (13)

$${\sigma }^{\left(6\right)}\left(x\right)-{\left(k\alpha H\right)}^2{\sigma }^{\left(4\right)}\left(x\right)=0$$

(62)

$$\begin{array}{cc}\begin{array}{c}{[1+{\left(k\alpha H\right)}^2\sigma ″\prime \left(x\right)-{\sigma }^{\left(5\right)}\left(x\right)]}_{x=\zeta }=0\hfill \\ \multicolumn{1}{c}{{[{\left(k\alpha H\right)}^2\sigma \prime \left(x\right)-\sigma ″\prime \left(x\right)]}_{x=\zeta }=0}\\ \multicolumn{1}{c}{{\sigma \prime \left(x\right)|}_{x=\zeta }=0}\end{array}\hfill & \begin{array}{c}{[-{\left(k\alpha H\right)}^2\sigma ″\left(x\right)+{\sigma }^{\left(4\right)}\left(x\right)]}_{x=\zeta }=0\hfill \\ \multicolumn{1}{c}{{[-{\left(k\alpha H\right)}^2\sigma \left(x\right)+\sigma ″\left(x\right)]}_{x=\zeta }=0}\\ \multicolumn{1}{c}{{\sigma \left(x\right)|}_{x=\zeta }=0}\end{array}\hfill \end{array}$$

(63)

which yields

$$\sigma (x,\zeta )=\frac{\sinh \left[k\alpha H(x-\zeta )\right]}{{\left(k\alpha H\right)}^5}-\frac{x-\zeta }{{\left(k\alpha H\right)}^4}-\frac{{(x-\zeta )}^3}{6{\left(k\alpha H\right)}^2}$$

(64)

Considering $L[w(\zeta )]=0$ , we can obtain the initial solution

$$\begin{array}{l}{w}_0\left(\zeta \right)=C_1+C_2\zeta +C_3{\zeta }^2+C_4{\zeta }^3\\ +C_5\cosh \left(k\alpha H·\zeta \right)+C_6\sinh \left(k\alpha H·\zeta \right)\end{array}$$

(65)

Because $N[w(t)]=0$ , the accurate analytical solution can be determined according to equation (25) (the starting point of the integral is chosen as $\zeta =1$ )

$$w\left(\zeta \right)=w_0\left(\zeta \right)+\underset{1}{\overset{\zeta }{\int }}\sigma (x,\zeta )·G\left(\zeta \right)\mathrm{d}x$$

(66)

where the expressions for $\sigma (x,\zeta )$ and $G(\zeta )$ are given in equations (64) and (61), respectively.

A linear system of equations with respect to the undetermined coefficients $C_1,C_2,\dots ,C_6$ can be obtained by substituting equation (66) into the boundary conditions in equation (59)

$$\left[\begin{array}{cccccc}0& 0& \frac{2\mathrm{sech}k\alpha H}{{\left(k\alpha H\right)}^2}& \frac{6\mathrm{sech}k\alpha H}{{\left(k\alpha H\right)}^2}& 1& \tanh k\alpha H\\ 0& 0& 0& 0& 1& \tanh k\alpha H\\ 0& 0& 0& -1& 0& 0\\ 1& 0& 0& 0& 1& 0\\ 0& -1& \frac{2}{{\left(\beta H\right)}^2}& 0& \frac{{\left(k\alpha H\right)}^2}{{\left(\beta H\right)}^2}& -k\alpha H\\ 0& 0& -\frac{2{\left(k\alpha H\right)}^2}{k^2{\left(\beta H\right)}^2}& 6& -\frac{{\left(k\alpha H\right)}^4}{k^2{\left(\beta H\right)}^2}& {\left(k\alpha H\right)}^3\end{array}\right]\left[\begin{array}{c}{C}_1\\ C_2\\ C_3\\ C_4\\ C_5\\ C_6\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\\ b_4\\ b_5\\ b_6\end{array}\right]$$

(67)

where

$$\begin{array}{c}{b}_1=0\hfill \\ b_2=\frac{H^4\mathrm{sech}k\alpha H}{{\left(k\alpha H\right)}^4\mathrm{EI}}·{p\left(\zeta \right)|}_{\zeta =1}\hfill \\ b_3=\frac{H^4}{6{\left(k\alpha H\right)}^2\mathrm{EI}}[{\frac{\mathrm{d}p\left(\zeta \right)}{\mathrm{d}\zeta }|}_{\zeta =1}+{\left(k\alpha H\right)}^2(1-k^{-2})\frac{P}{H}]\hfill \\ b_4=-\underset{1}{\overset{0}{\int }}\sigma (x,0)·G\left(x\right)\mathrm{d}x\hfill \\ b_5=-\underset{1}{\overset{0}{\int }}[\frac{1}{{\left(\beta H\right)}^2}·{\frac{{\mathrm{d}}^2\sigma (x,\zeta )}{\mathrm{d}{\zeta }^2}|}_{\zeta =0}-{\frac{\mathrm{d}\sigma (x,\zeta )}{\mathrm{d}\zeta }|}_{\zeta =0}]·G\left(x\right)\mathrm{d}x\hfill \\ b_6=-\frac{H^4}{\mathrm{EI}}\underset{0}{\overset{1}{\int }}p\left(\zeta \right)\mathrm{d}\zeta \hfill \\ -\underset{1}{\overset{0}{\int }}[{\frac{{\mathrm{d}}^3\sigma (x,\zeta )}{\mathrm{d}{\zeta }^3}|}_{\zeta =0}-\frac{{\left(k\alpha H\right)}^2}{k^2{\left(\beta H\right)}^2}{·\frac{{\mathrm{d}}^2\sigma (x,\zeta )}{\mathrm{d}{\zeta }^2}|}_{\zeta =0}]·G\left(x\right)\mathrm{d}x\hfill \end{array}$$

(68)

The values of $C_1,C_2,\dots ,C_6$ can be acquired by solving the linear system represented by equation (67) and the solution of the differential equation (equation (56)) can be obtained by substituting the values of $C_1,C_2,\dots ,C_6$ in equations (65) and (66). As mentioned in the previous section, the integrals in equation (66) will be calculated herein by a numerical approach (numerical integration using the trapezoidal rule or Simpson’s rule) for the sake of practicality and convenience. To avoid a numerical differential, the differential can be calculated according to the property of the generalized Lagrange multiplier (equation (17)) as follows

$$w^{\left(m\right)}\left(\zeta \right)=w_0^{\left(m\right)}\left(\zeta \right)+\underset{1}{\overset{\zeta }{\int }}\frac{{\partial }^{\left(m\right)}\sigma (x,\zeta )}{\partial x^{\left(m\right)}}·q\left(\zeta \right)\mathrm{d}x$$

(69)

A computer program based on the above derivation named ODETB (solver for an ODE of a tall building structure under static loads) was developed using Julia Language, ³¹ which is a high-level, high-performance, dynamic programming language for technical computing. This ODETB is capable of solving ODEs (such as equation (56)) for tall building structures under arbitrary lateral loads.

To verify the ODETB, a simple example is presented wherein a tall building is carrying a distributed inverted triangle load $p(\zeta )=\zeta ·\mathrm{EI}/H^4$ , thus making it easy to obtain an analytical solution. The key parameters of the tall building are set as $k^2=1.2$ and $k\alpha H=\beta H=10.0$ . A comparison between the results obtained by the semi-analytical solution (ODETB) and analytical solution is shown in Table 2, in which the relative error is calculated by

$$\mathrm{E}=\frac{\underset{0}{\overset{1}{\int }}|\tilde{w}\left(\zeta \right)-w\left(\zeta \right)|\mathrm{d}\zeta }{\underset{0}{\overset{1}{\int }}\left|w\left(\zeta \right)\right|\mathrm{d}\zeta }$$

(70)

where $w(\zeta )$ represents the analytical solution and $\tilde{w}(\zeta )$ represents the numerical solution. The error data listed in Table 2 indicate that ODETB exhibits relatively high accuracy even if there are only few points in the solution domain for numerical integration, and the suggested number of points used in the numerical integration should be 50 considering the computational efficiency and accuracy.

Table 2.

Relative errors of ODETB results compared with analytical solutions.

Number of points used in numerical integration	Relative error
	Displacement	Rotation angle	Curvature
3	6.90E−02	3.00E−01	2.80E+00
5	1.10E−03	3.30E−03	6.30E−02
9	4.70E−05	6.30E−05	2.90E−03
17	2.70E−06	1.60E−06	1.70E−04
33	1.60E−07	1.20E−07	1.00E−05
65	9.90E−09	9.30E−09	6.10E−07
129	4.00E−10	3.90E−10	3.80E−08
257	1.90E−10	2.30E−10	2.30E−09
513	2.30E−10	2.90E−10	4.20E−10
1025	2.30E−10	2.80E−10	4.70E−10

Solution of sixth-order ODE eigenvalue problem

Comparing equation (58) with equations (10) and (29), the operators can be identified as

$$\begin{array}{c}\begin{array}{cc}L\left[y\left(\zeta \right)\right]=\frac{{\mathrm{d}}^{6}y}{\mathrm{d}{\zeta }^6}-{\left(k\alpha H\right)}^2\frac{{\mathrm{d}}^{4}y}{\mathrm{d}{\zeta }^4}\hfill & G\left(\zeta \right)=0\hfill \end{array}\hfill \\ N\left[y\left(\zeta \right)\right]=-{\left(\lambda H\right)}^4[\frac{{\mathrm{d}}^{2}y}{\mathrm{d}{\zeta }^2}-{\left(k\alpha H\right)}^2(1-k^{-2})y]\hfill \\ \overline{N}\left[y\left(\zeta \right)\right]=-[\frac{{\mathrm{d}}^{2}y}{\mathrm{d}{\zeta }^2}-{\left(k\alpha H\right)}^2(1-k^{-2})y]\hfill \end{array}$$

(71)

The generalized Lagrange multiplier is same as that in equation (64), and the initial solution is also same as that in equation (65) because the operator L in the above equation is identical to the one in equation (61). To make the solution process simpler, the new initial solution is determined using the boundary conditions at $\zeta =1$ , thereby eliminating three of the undetermined coefficients

$$y_0\left(\zeta \right)=C_1{f}_{1,0}\left(\zeta \right)+C_2{f}_{2,0}\left(\zeta \right)+C_3{f}_{3,0}\left(\zeta \right)$$

(72)

where

$$\begin{array}{c}{f}_{1,0}\left(\zeta \right)=\sinh \left[k\alpha H(\zeta -1)\right]-k\alpha H(\zeta -1)\hfill \\ \begin{array}{c}{f}_{2,0}\left(\zeta \right)={\left(\lambda H\right)}^4\hfill \\ \multicolumn{1}{c}{[\cosh \left[k\alpha H(\zeta -1)\right]-0.5{\left(k\alpha H\right)}^2{(\zeta -1)}^2-1]+{\left(k\alpha H\right)}^4}\\ \multicolumn{1}{c}{f_{3,0}\left(\zeta \right)={\left(\lambda H\right)}^4{(\zeta -1)}^3-6{\left(k\alpha H\right)}^2(\zeta -1)}\end{array}\hfill \end{array}$$

(73)

Because there is an unknown variable $(\lambda H{)}^4$ in the above expressions, the initial solution should be rearranged as

$$\begin{array}{cc}\multicolumn{1}{c}{y_0\left(\zeta \right)=}& C_1{g}_{1,0}\left(\zeta \right)+C_2[{\left(\lambda H\right)}^4{g}_{2,0}\left(\zeta \right)+g_{3,0}\left(\zeta \right)]\hfill \\ \multicolumn{1}{c}{}& +C_3[{\left(\lambda H\right)}^4{g}_{4,0}\left(\zeta \right)+g_{5,0}\left(\zeta \right)]\hfill \end{array}$$

(74)

where

$$\begin{array}{c}{g}_{1,0}\left(\zeta \right)=\sinh \left[k\alpha H(\zeta -1)\right]-k\alpha H(\zeta -1)\hfill \\ \multicolumn{1}{c}{\begin{array}{c}{g}_{2,0}\left(\zeta \right)=[\cosh \left[k\alpha H(\zeta -1)\right]-0.5{\left(k\alpha H\right)}^2{(\zeta -1)}^2-1]\hfill \\ \multicolumn{1}{c}{g_{3,0}\left(\zeta \right)={\left(k\alpha H\right)}^4}\\ \multicolumn{1}{c}{g_{4,0}\left(\zeta \right)={(\zeta -1)}^3}\\ \multicolumn{1}{c}{g_{5,0}\left(\zeta \right)=-6{\left(k\alpha H\right)}^2(\zeta -1)}\end{array}}\end{array}$$

(75)

The ith iterative solution can then be determined on the basis of equations (30) and (74) as

$$\begin{array}{cc}\multicolumn{1}{c}{y_i\left(\zeta \right)=}& C_1\left[\sum _{l=0}^i{\left(\lambda H\right)}^{4l}{g}_{1,l}\left(\zeta \right)\right]\hfill \\ \multicolumn{1}{c}{}& +C_2\sum _{l=0}^i{\left(\lambda H\right)}^{4l}[{\left(\lambda H\right)}^4{g}_{2,l}\left(\zeta \right)+g_{3,l}\left(\zeta \right)]\hfill \\ \multicolumn{1}{c}{}& +C_3\sum _{l=0}^i{\left(\lambda H\right)}^{4l}[{\left(\lambda H\right)}^4{g}_{4,l}\left(\zeta \right)+g_{5,l}\left(\zeta \right)]\hfill \end{array}$$

(76)

where

$$g_{k,l}\left(\zeta \right)=-{\displaystyle \underset{1}{\overset{\zeta }{\int }}\sigma \left(x,\zeta \right)·\left[{g^{″}}_{k,l-1}\left(x\right)-{\left(k\alpha H\right)}^2\left(1-k^{-2}\right)g_{k,l-1}\left(x\right)\right]\mathrm{d}x}$$

(77)

The differential involved in equations (77) and (60) can be computed as follows according to equation (17)

$$g_{k,l}^{\left(m\right)}\left(\zeta \right)=-{\displaystyle \underset{1}{\overset{\zeta }{\int }}\frac{{\partial }^m\sigma \left(x,\zeta \right)}{\partial x^m}·\left[{g^{″}}_{k,l-1}\left(x\right)-{\left(k\alpha H\right)}^2\left(1-k^{-2}\right)g_{k,l-1}\left(x\right)\right]\mathrm{d}x}$$

(78)

where m is the order of the differential and m < 6.

A linear system of equations with respect to the undetermined coefficients $C_1,C_2,C_3$ can be obtained by substituting the boundary value calculated by equation (76) into the boundary conditions in equation (60)

$$\left[\begin{array}{ccc}{p}_{11}\hfill & p_{12}\hfill & p_{13}\hfill \\ \multicolumn{1}{c}{p_{21}}& p_{22}\hfill & p_{23}\hfill \\ \multicolumn{1}{c}{p_{31}}& p_{32}\hfill & p_{33}\hfill \end{array}\right]\left\{\begin{array}{c}{C}_1\hfill \\ \multicolumn{1}{c}{C_2}\\ \multicolumn{1}{c}{C_3}\end{array}\right\}=\left\{\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\right\}$$

(79)

where

$$\begin{array}{c}{p}_{11}=\sum _{i=0}^n{\left(\lambda H\right)}^{4i}{g}_{1,i}\left(0\right)\hfill \\ p_{12}=\sum _{i=0}^n{\left(\lambda H\right)}^{4i}[{\left(\lambda H\right)}^4{g}_{2,i}\left(0\right)+g_{3,i}\left(0\right)]\hfill \\ p_{13}=\sum _{i=0}^n{\left(\lambda H\right)}^{4i}[{\left(\lambda H\right)}^4{g}_{4,i}\left(0\right)+g_{5,i}\left(0\right)]\hfill \\ p_{21}=\sum _{i=0}^n{\left(\lambda H\right)}^{4i}[\frac{1}{{\left(\beta H\right)}^2}·{g″}_{1,i}\left(0\right)-{g\prime }_{1,i}\left(0\right)]\hfill \\ p_{22}=\sum _{i=0}^n{\left(\lambda H\right)}^{4i}[\frac{1}{{\left(\beta H\right)}^2}·[{\left(\lambda H\right)}^4{g″}_{2,i}\left(0\right)+{g″}_{3,i}\left(0\right)]-[{\left(\lambda H\right)}^4{g\prime }_{2,i}\left(0\right)+{g\prime }_{3,i}\left(0\right)]]\hfill \\ p_{23}=\sum _{i=0}^n{\left(\lambda H\right)}^{4i}[\frac{1}{{\left(\beta H\right)}^2}·[{\left(\lambda H\right)}^4{g″}_{4,i}\left(0\right)+{g″}_{5,i}\left(0\right)]-[{\left(\lambda H\right)}^4{g\prime }_{4,i}\left(0\right)+{g\prime }_{5,i}\left(0\right)]]\hfill \end{array}$$

(80)

$$\begin{array}{c}{p}_{31}=\sum _{i=0}^n{\left(\lambda H\right)}^{4i}[g_{1,i}″\prime \left(0\right)-\frac{{\left(k\alpha H\right)}^2}{k^2}·g_{1,i}\prime \left(0\right)+{\left(\lambda H\right)}^4\underset{0}{\overset{1}{\int }}{g}_{1,i}\left(\zeta \right)\mathrm{d}\zeta ]\hfill \\ p_{32}=\sum _{i=0}^n{\left(\lambda H\right)}^{4i}\left[\begin{array}{c}{\left(\lambda H\right)}^8\underset{0}{\overset{1}{\int }}{g}_{2,i}\left(\zeta \right)\mathrm{d}\zeta +{\left(\lambda H\right)}^4[g_{2,i}″\prime \left(0\right)-\frac{{\left(k\alpha H\right)}^2}{k^2}·g_{2,i}\prime \left(0\right)]\hfill \\ \multicolumn{1}{c}{+{\left(\lambda H\right)}^4\underset{0}{\overset{1}{\int }}{g}_{3,i}\left(\zeta \right)\mathrm{d}\zeta +g_{3,i}″\prime \left(0\right)-\frac{{\left(k\alpha H\right)}^2}{k^2}·g_{3,i}\prime \left(0\right)}\end{array}\right]\hfill \\ p_{33}=\sum _{i=0}^n{\left(\lambda H\right)}^{4i}\left[\begin{array}{c}{\left(\lambda H\right)}^8\underset{0}{\overset{1}{\int }}{g}_{4,i}\left(\zeta \right)\mathrm{d}\zeta +{\left(\lambda H\right)}^4[g_{4,i}″\prime \left(0\right)-\frac{{\left(k\alpha H\right)}^2}{k^2}·g_{4,i}\prime \left(0\right)]\hfill \\ \multicolumn{1}{c}{+{\left(\lambda H\right)}^4\underset{0}{\overset{1}{\int }}{g}_{5,i}\left(\zeta \right)\mathrm{d}\zeta +g_{5,i}″\prime \left(0\right)-\frac{{\left(k\alpha H\right)}^2}{k^2}·g_{5,i}\prime \left(0\right)}\end{array}\right]\hfill \end{array}$$

(81)

According to Cramer’s rule, the following equation must hold true if the above linear system has non-zero solutions

$$\left|\begin{array}{ccc}{p}_{11}\hfill & p_{12}\hfill & p_{13}\hfill \\ \multicolumn{1}{c}{p_{21}}& p_{22}\hfill & p_{23}\hfill \\ \multicolumn{1}{c}{p_{31}}& p_{32}\hfill & p_{33}\hfill \end{array}\right|=0$$

(82)

It can be predicted that equation (82) can be arranged as a polynomial equation of $(\lambda H{)}^4$ of order 3(i + 1), where i is the iteration count. The positive real roots of the equation are simply the eigenvalues of the ODE given by equation (58). The eigenmodes can then be abstracted by substituting the eigenvalues into equations (74) and (79), and then eliminating two of the three undetermined coefficients $C_1,C_2,\mathrm{and}{C}_3$ .

Based on the above discussions, we have developed a corresponding computer program using Julia to solve the ODE eigenvalue problems of tall buildings by the semi-analytical VIM. To validate this program, which is called EIGTB (solver for a ODE eigenvalue problem of a tall building structure for free vibration), two examples referred to by Chai and Chen ³² were recomputed using EIGTB and two other approximate methods: the equivalent sandwich beam method ³³ and Rutenberg’s ³⁴ method. The numerical results of the eigenvalues of these examples are listed in Table 3, and the eigenmodes are shown in Figure 1 (EIGTB results were obtained using 50 integral points, involving only 6 iterations). These results confirm that the semi-analytical VIM and EIGTB program proposed in this article are valid and accurate.

Table 3.

Eigenvalues of continuum tall building structures obtained by various methods.

Eigenvalues $(\lambda H{)}^4$		First order	Second order	Third order	Fourth order
Example 1:
$\begin{array}{c}{k}^2=1.2205\hfill \\ k\alpha H=4.9505\hfill \\ \beta H\to \infty \hfill \end{array}$	Chai and Chen 32	43.699	871.80	5129.7	17218
	EIGTB	43.6987	871.8024	5129.6478	17422.4116
	Sandwich beam 33	38.6642	856.3319	4962.2431	16959.1977
	Rutenberg 34	43.4881	961.4456	5217.8186	17360.2280
Example 2:
$\begin{array}{c}{k}^2=1.2729\hfill \\ k\alpha H=17.95727\hfill \\ \beta H\to \infty \hfill \end{array}$	Chai and Chen 32	55.032	1746	10945	33939
	EIGTB	55.0318	1746.0435	10944.6228	33805.3503
	Sandwich beam 33	54.6012	1837.1746	11176.5327	34103.1317
	Rutenberg 34	54.8940	1873.2208	11559.1689	35427.6689

Figure 1.

Eigenmodes of continuum tall building structures solved by EIGTB: (a) example 1 and (b) example 2.

Conclusion

In this study, the original VIM is simplified and modified into a semi-analytical method to solve high-order complicated ODEs in a more convenient and practical manner.

1. The general form of the differential equation with respect to the general Lagrange multiplier is explicitly stated, which facilitates the identification of the general Lagrange multiplier.

2. For linear ODEs, the operator L can be eliminated during the iterative solving procedure, which reduces the iterative computational burden.

3. The integrals involved in the iterations are suggested to be calculated numerically for practical use.

4. The simplified VIM has been proven to be an efficient and practical method for solving complicated ODEs (including eigenvalue ODEs), and the implementation of this method is both simple and convenient.