In this thesis, the analysis of the method of fundamental solution(MFS) is expanded for biharmonic equations. The bounds of errors are derived for the traditional and the Almansi's approaches in bounded simply-connected domains. The exponential and the polynomial convergence rates are obtained from highly and finite smooth solutions, respectively. Also the bounds of condition number are derived for the disk domains, to show the exponential growth rates. The analysis in this thesis is the first time to provide the rigor analysis of the CTM for biharmonic equations, and the intrinsic nature of accuracy and stability is similar to that of Laplace's equation. Numerical experiment are carried out for both smooth and singularity problems. The numerical results coincide with the theoretical analysis made. When the particular solutions satisfying the biharmonic equation can be found, the method of particular solutions(MPS) is always superior to MFS, supported by numerical examples. However, if such singular particular solutions near the singular points can not be found, the local refinement of collocation nodes and the greedy adaptive techniques can be used. It seems that the greedy adaptive techniques may provide a better solution for singularity problems. Beside, the numerical solutions by Almansi's approaches are slightly better in accuracy and stability than those by the traditional FS. Hence, the MFS with Almansi's approaches is recommended, due to the simple analysis, which can be obtained directly from the analysis of MFS for Laplace's equation.