## 3 The Model's Difference Equation and Its Solution  
To solve the above model, the pile shaft is divided into \(n\) units vertically, each of length \(\lambda\). For convenience, four virtual units are added at the top and bottom of the pile, and the unit and node numbers are shown in Figure 1. For any pile node \(i\), the difference equation of node \(i\) can be obtained by substituting a one-dimensional central difference formula into Equation (1), according to the established pile deflection differential equation. \[y_{i+2}-\left(4-\frac{N\lambda^{2}}{EI}\right)\cdot y_{i+1}+\left(6-\frac{2N \lambda^{2}}{EI}+\frac{b_{0}c_{zi}\lambda^{4}}{EI}\right)\cdot y_{i}-\left(4- \frac{N\lambda^{2}}{EI}\right)\cdot y_{i-1}+y_{i-2}-\frac{q_{i}\lambda^{4}d}{ EI}=0 \tag{19}\]  
\[\text{Substituting }4-\frac{N\lambda^{2}}{EI}=k,\frac{b_{0}c_{zi}\lambda^{4}}{ EI}=A_{i}\text{ and }\frac{\lambda^{4}d}{EI}=G\text{ into Equation (\ref{eq:pile_def}), Equation (\ref{eq:pile_def}) can be written as  
\[y_{i+2}-k\cdot y_{i+1}+(2k-2+A_{i})\cdot y_{i}-k\cdot y_{i-1}+y_{i-2}-G\cdot q _{i}=0 \tag{20}\]  
Corresponding rotation (\(\theta_{i}\)), bending moment (\(M_{i}\)), and shear force (\(Q_{i}\)) can be obtained as  
\[\theta_{i}=\frac{1}{2\lambda}\cdot(y_{i+1}-y_{i-1}) \tag{21}\]  
\[M_{i}=\frac{EI}{\lambda^{2}}\cdot(y_{i+1}-2y_{i}+y_{i-1}) \tag{22}\]  
\[Q_{i} =\frac{dM}{dz}+N\cdot\frac{dy}{dz} \tag{23}\] \[=\frac{EI}{2\lambda^{3}}\cdot(y_{i+2}-2y_{i+1}+2y_{i-1}-y_{i-2})+ \frac{N}{2\lambda}\cdot(y_{i+1}-y_{i-1})\]  
An equilibrium equation such as Equation (19) can be obtained for any pile shaft node \(i\). There are \(n+1\) equations. To solve the system, four supplementary equations need to be determined from the boundary conditions. Two equations can be derived using the known active load of the pile top:  
\[M_{3}=\frac{EI}{\lambda^{2}}\cdot(y_{4}-2y_{3}+y_{2})=M_{0} \tag{24}\]  
\[Q_{3}=\frac{EI}{2\lambda^{3}}\cdot(y_{5}-2y_{4}+2y_{2}-y_{1})+\frac{N}{2 \lambda}\cdot(y_{4}-y_{2})=H_{0} \tag{25}\]  
and the other two can be derived according to the boundary conditions for the pile bottom. When the pile bottom is not embedded in the rock stratum, it can be regarded as free. That is, the shear force at node \(n+3\) is 0. When the pile bottom is supported in the soil layer and \(al\geq 2.5\) or the pile end is supported in the rock stratum and \(al\geq 3.5\), the section deformation of the pile bottom is very small. Thus, the bending moment at the pile bottom can treated as 0:  
\[M_{n+3}=\frac{EI}{\lambda^{2}}\cdot(y_{n+4}-2y_{n+3}+y_{n+2})=0 \tag{26}\]  
\[y_{i+2}-\left(4-\frac{N\lambda^{2}}{EI}\right)\cdot y_{i+1}+\left(6-\frac{2N \lambda^{2}}{EI}+\frac{b_{0}c_{zi}\lambda^{4}}{EI}\right)\cdot y_{i}-\left(4- \frac{N\lambda^{2}}{EI}\right)\cdot y_{i-1}+y_{i-2}-\frac{q_{i}\lambda^{4}d}{ EI}=0 \tag{19}\]0When the pile end is embedded in the rock stratum, it is regarded as fixed, and in this case, the section angle and displacement at node \(n+3\) are 0:  
\[y_{i+2}-\left(4-\frac{N\lambda^{2}}{EI}\right)\cdot y_{i+1}+\left(6-\frac{2N \lambda^{2}}{EI}+\frac{b_{0}c_{zi}\lambda^{4}}{EI}\right)\cdot y_{i}-\left(4- \frac{N\lambda^{2}}{EI}\right)\cdot y_{i-1}+y_{i-2}-\frac{q_{i}\lambda^{4}d}{ EI}=0 \tag{19}\]1  
\[y_{i+2}-\left(4-\frac{N\lambda^{2}}{EI}\right)\cdot y_{i+1}+\left(6-\frac{2N \lambda^{2}}{EI}+\frac{b_{0}c_{zi}\lambda^{4}}{EI}\right)\cdot y_{i}-\left(4- \frac{N\lambda^{2}}{EI}\right)\cdot y_{i-1}+y_{i-2}-\frac{q_{i}\lambda^{4}d}{ EI}=0 \tag{19}\]2  
Here, the equations are formed by combining the different pile bottom boundary conditions, pile-top load conditions, and pile node difference equations. After solving this system of equations, all node displacements \(y_{i}\) can be obtained successively from \(y_{0}\). Then, the rotation, bending moment, and shear force of each node of the pile can be obtained by Equation (21) through Equation (23). In fact, since the horizontal displacement (\(y_{i}\)) of pile at point \(i\) is unknown, \(c_{zi}\) cannot be determined from the \(p\)-\(y\) curve.