If you are working with finite difference methods, Any equation PDE can be written as Pu=f where P is an operator. We introduce P_{k,h} the approximation operator for u and R_{k,h} the approximation operator of f, then any finite difference scheme can be written P_{k,h} v= R_{k,h} f . Now if you take any smooth function $Phi(t,x)$ the order of accuracy is (p,q) where P_{k,h}(Phi)-R_{k,h}(P(Phi))=O(k^p+h^q)
So take any smmoth function and compute (P_{k,h}(Phi)-R_{k,h}(P(Phi)))/(k^p+h^q) if this quantity is constant or less than a constant then the order is (p,q)
Roughly, one may get idea through halving the step length. Suppose N1 is a solution for step length h, and N2 is another solution at h/2. Then N1 + O(hp) = N2 + O (hp/ 2p). Check how | N1 - N2| = hp (1- 1/ 2 p) is varying.
It is of interest to see that there are error estimates in a cumulative per step estimation but that there are errors in introducing parasitic roots that differ markedly from the solution to the original PDE. This may show up as an instability that masks the true solution by creating a system response from the integration methodology.
There is a formal method for carrying out such analysis using so-called modified equations to determine the pseudo-physical side effects that various differencing schemes introduce into a numerical solution (e.g. how much "artificial dissipation" arises). Tannehill, Anderson, and Pletcher "Computational fluid mechanics and heat transfer" (1997) has a very thorough discussion of the method. This would be the general approach.
The "method of manufactured solutions" mentioned by Aldo above is also good (though less general) because it furnishes a direct basis of numerical comparison between an exact and a numerical solution (where the actual physical problem is not terribly important). As a shameless plug, Wendl (1999) "A method for obtaining benchmark Navier–Stokes solutions" Int. J. Num. Meth. Fluids 31(3), 657-660 describes a general method for building fairly sophisticated solutions of this type.
If you used analytic-approximate methods for examples:ADM, VIM, HPM and HAM...etc, you can simple calculate the residual error(i.e. error remainder). However, for numerical methods, you can solve your problem by say RK4 OR Euler methods using bulit in functions in MATLAB, MATHEMATICA and considering it as exact, then evaluate the relative errors or RMS, some of above answers are also very useful.
You have someways to do that when you do not have access to exact solution:
1- Make a exact solution (it is called method of manufactured solution (MMS) or Prescribed Solution Forcing (PSF)). That is not easy and fun but works very well!
2- Do cross-code verification: Solve the problem on the pre verified solver with high resolution and then use that to compare your results and find formal order of convergence.
3- Compare versus dense mesh solution of the same numerical solver! that usually hits practical limitation in 2D and 3D but in theory you can do that.
4- Do Richardson Extrapolation to figure out the order of convergence. It has some practical issue and consideration but this is the most common answer to your question.
A) I wrote a python subroutine which can does the above for you, see here
Vavuq.org
B) See chapter 5 of this report for examples and details of implementation:
On the issues associated with 2-D modeling for flood mapping purposes
First thing is if constraints to the PDE equations are unknown then it is difficult to assume any function that will satisfy the governing equilibrium equation and constraints to the equation. There you can assume classical methods like Galerkin,Rayleigh-Ritz,etc