In this thesis, we consider the numerical approximation of high order geometric Partial Differential Equations (PDEs). We first consider high order PDEs defined on surfaces in the 3D space that are represented by single-patch tensor product NURBS. Then, we spatially discretize the PDEs by means of NURBS-based Isogeometric Analysis (IGA) in the framework of the Galerkin method. With this aim, we consider the construction of periodic NURBS function spaces with high degree of global continuity, even on closed surfaces. As benchmark problems for the proposed discretization, we propose Laplace-Beltrami problems of the fourth and sixth orders, as well as the corresponding eigenvalue problems, and we analyze the impact of the continuity of the basis functions on the accuracy as well as on computational costs. The numerical solution of two high order phase field problems on both open and closed surfaces is also considered: the fourth order Cahn-Hilliard equation and the sixth order crystal equation, both discretized in time with the generalized-alpha method. We then consider the numerical approximation of geometric PDEs, derived, in particular, from the minimization of shape energy functionals by L^2-gradient flows. We analyze the mean curvature and the Willmore gradient flows, leading to second and fourth order PDEs, respectively. These nonlinear geometric PDEs are discretized in time with Backward Differentiation Formulas (BDF), with a semi-implicit formulation based on an extrapolation of the geometry, leading to a linear problem to be solved at each time step. Results about the numerical approximation of the two geometric flows on several geometries are analyzed. Then, we study how the proposed mathematical framework can be employed to numerically approximate the equilibrium shapes of lipid bilayer biomembranes, or vesicles, governed by the Canham-Helfrich curvature model. We propose two numerical schemes for enforcing the conservation of the area and volume of the vesicles, and report results on benchmark problems. Then, the approximation of the equilibrium shapes of biomembranes with different values of reduced volume is presented. Finally, we consider the dynamics of a vesicle, e.g. a red blood cell, immersed in a fluid, e.g. the plasma. In particular, we couple the curvature-driven model for the lipid membrane with the incompressible Navier-Stokes equations governing the fluid. We consider a segregated approach, with a formulation based on the Resistive Immersed Surface method applied to NURBS geometries. After analyzing benchmark fluid simulations with immersed NURBS objects, we report numerical results for the investigation of the dynamics of a vesicle under different flow conditions.