High energy astrophysics is a relatively new discipline that has witnessed tremendous advancement in the past 25 years or so with the advent of space technology, which provides the field with an indispensable observational tool in the X-ray and γ-ray range previously unreachable by ground-based observations. Ongoing discoveries in this field demand new theories to explain new phenomena. In this thesis, I summarize my work on one of these topics in high energy astrophysics, the theory of accretion disks around black holes. Accretion disk models are applied to the dynamics of compact objects surrounded by material with organized angular momentum. The disk forms due to the process by which differentially rotating fluid layers interact with each other and transport angular momentum outwards with the gas spiraling into the central object, a black hole in this case. This process, known as accretion, is a very effective way to tap the energy of the gravitational field. In fact, for a non-rotating black hole, the conversion rate from mass to energy is about six percent. This is about ten times as effective as the sun, while for an extreme rotating black hole, as much as forty percent of the rest mass of the accreted gas is converted into energy. Based on the X-ray spectra and time variability of black hole candidates, I have developed a so-called radial-zoned thermal accretion disk model by dividing the disk radially into three zones. The outer zone is optically thick and the inner and middle zones are optically thin, with the latter being irradiated by the outer part of the disk. This model can produce a composite spectrum of black body radiation from the outer zone, optically thin comptonized thermal bremsstrahlung radiation from the inner zone, and inverse comptonized soft photon radiation from the middle. In this work, I have studied both the linear and nonlinear dynamic evolution and stability of this accretion disk model by solving the whole set of partial differential equations numerically. Spectral fitting and time series analysis to high energy satellite data have been done based on the steady disk model to fix parameters and to cross-examine the results from other methods.