I have discrete data of a 2D function defined as

theta = linspace(0,pi,nTheta);
phi = linspace(0,2*pi,nPhi);
p=zeros(nPhi,nTheta);%only to show the dimension of my matrix
[np,nt]=ndgrid(phi,theta);
f1 = griddedInterpolant(np,nt,p,'spline');
f2= @(np,nt) f1(np,nt);
integral2(f2,0,2*pi,0,pi)

Note that p is calculated from a complex physical problem, but i showed above how it is initialized.

Also, I can increase nTheta and nPhi, which leads to more accurate calculation of p.

My calculated function (with nPhi=400,nTheta=200) is something like:

I tried 3 ways :

• using Trapz function
• using the code above but with linear interpolation for gridded interpolant
• using the code above with spline interpolation

Although the spline is better than others, i still need to increase nPhi and nTheta, which makes it impossible for me to do the simulation due to its cost. Is there any suggestion except these 3 methods or any general suggestion how i can do this computation more efficient? (I also took advantage of the symmetry in both directions)

Note that the shape of my function varies in each time step, so a local mesh refinement might be challenging because i don't know the detail of my function in advance.