This dissertation investigates the approximation techniques and instability for nonlinear one-link and two-link flexible manipulators operating in a vertical plane. The flexible components of the arms are modeled by Euler-Bernoulli beam theory, and the nonlinearity arises due to gravitation and large angle rotation. Two different methods are used for describing the deformation of the flexible members. They are the infinite-dimensional coordinates and finite-dimensional coordinates. Variational principle (Hamilton's principle) is used to generate the governing differential equations and boundary conditions. The distributed coordinate approach results in partial differential equations and the Ritz approximation leads to nonlinear ordinary differential equations. The two sets of differential equations are developed in parallel to give a comparison of the static solution and frequency domain characteristics. Furthermore, the full nonlinear ordinary differential equations is integrated forward numerically for the dynamic response to maneuver transient reference inputs;Among many interesting topics, the following issues are studied (1) Use of the exact solution as the benchmark for the approximation solutions. (2) Effect of the admissible comparison functions to accuracy of approximation solution. (3) Influence of set point selection to the linearized open- and closed-loop flexible dynamic systems. (4) Effect of linear visco-elastic damping to the Laplace transform domain and time domain behavior of a flexible arm. (5) Effects of gravitational force on the closed-loop control systems with PID control at each joint;Various comparison functions are used to discretize the equations of motion of the deformable arm. Since the Ritz method requires only the essential boundary conditions to be satisfied and places no restriction on the natural boundary conditions, it allows the use of many different types of shape functions. Among the various sets of possible shape functions, only certain sets would satisfy both essential and natural boundary conditions, while the rest satisfy only the geometrical boundary conditions. Examples are given to show the importance of selecting comparison functions. In these examples both the exact and approximate solutions are obtained either in closed form or numerically. The effect of the discretization is analyzed in the Laplace transform domain by comparing the approximate solutions with the closed-form solution, and the causes of the differences in results are identified and analyzed. Finally, Liapunov's direct method is used to re-examine the stability characteristics. A sufficient condition for a stable PD control system is derived.