Evaluating integrals multiple input functions and a matrix of data

Dear all;

I'm working on a project and have need to numerically calculate an overlap integral over r,theta, phi. The function I'm looking at however is also a spherical harmonic (of specific form with l,m). Each coordinate input is a 3D meshgrid as this is being evaluated over a volume.

Currently my integrated function is of the form;

function rY_nm = (R.^n).*spharm(n,m,Theta,Phi); <-- 25ish lines + plotting.
assign variable: rY_nm(n,m,Theta,Phi);
funintegrall = funtest .*rY_nm.*sin(Theta);
inttheta2 = trapz(theta,fint); intphi2 = trapz(phi,inttheta2); intr1(n,mi) = trapz(r,intphi2);

Within a for loop which runs over n,m resulting in a set of coefficients at a each specific intr1(n,m). Here r, theta, phi are the linear coordinates (i.e unmeshed setup).

For my work however I need to use a much higher accuracy numerical method (even at 500x500x500 points the answers aren't following the pattern expected). I have access to a few such as integral and d001ah from the NAG libaries, however these require a function handle as an input which I have zero experience with. I have spent a significant amount of time trying to assign a function handle e.g fint = (R, Theta, Phi) (R.^n) .*spharm(n,m,Theta,Phi) or fint = (R, Theta, Phi,n,m) (R.^n).*spharm(n,m,Theta,Phi) to little success. Would anyone be able to advise or direct me into how i would assign function handles in this type of system?

Answer 1

[Star Strider]

If you have a speific problem, we can likely provide much more specific solutions.

For the time being, see the documentation on Function Basics (link), and Anonymous Functions (link).

Function handles are not difficult to define and use in practice. If you have yoiur own function file, for example:

function y = myFunction(x,n)
    y = x.^n;
end

saved in your MATLAB path as ‘myFunction.m’, you would refer to it (for example in the integral (link) function) as:

n = 0:5;
a = 0;
b = 10;
int_y = integral(@(x)myFunction(x,n), a, b, 'ArrayValued',1);

since as an argument to another function, you would need to pass it as a function handle, using the ‘@’ operator.

This should at least help you understand the barest of essentials of function handle use, as well as an introduction to the integral function.

If you want significant accuracy in your integrals, I would use integral and related functions, if you have functions you want to integrate. The trapz function is for vectors and matrices.

[Logan Turton]

I'll try to be more specific here, apologies!

I have some set of matrix data representing a magnetic field (D) over a known spherical volume (i.e my coordinates R,Theta,Phi, each of which is a 3D matrix representing the coordinates over a volume).

Mathematically one can linearly decompose this into a set of spherical harmonics. The weightings of each coefficient present can then be found by performing the overlap integral of the function against a specific spherical harmonic in spherical coordinates.

As part of my testing and other work, I create an input function which should be mostly a single harmonic with others present giving me this 3D set of data. I have another function, [Y_nm] = spharm(n,m,Theta,Phi) which creates the spherical harmonics over the same coordinates, such rYnm = R.^n .*spharm(n,m,Theta,Phi).

I can create function handles for the R term and spharm thanks to your help, however my matrix D is still just a matrix.

Ynm = @(n,m, Theta,Phi) spharm(n,m,Theta,Phi);

Rn = @(R) R.^n;

rYnm = @(R,n,m,Theta,Phi) Rn.*Ynm;

D = size NxNxN

Trapezodial and simpsons integrals even with a large number of points are too inaccurate (main source: Linfield Numerical Methods Using MATLAB, plus some of my own testing) due to me working with nT fields.Hence the need to use something more accurate numerically such as the functionality provided in integral. It is however essential that I can still take the numerical volume of points to provide visualiations of my fields at any particular time in the calculations.

Any thoughts or suggestions on what I maybe able to do would be greatly appreciated.

My current approach looks something as follows:

Ynm = @(n,m,Theta,Phi) spharm(n,m,Theta,Phi);

Rn = @(R,n) R.^n;

rYnm = @(n,m,R,Theta,Phi) Rn(R,n).*Ynm(n,m,Theta,Phi);

funnm = @(R,Theta,Phi) fnm;

fint = @(n,m,R,Theta,Phi)...

funnm(R,Theta,Phi)*rYnm(n,m,R,Theta,Phi)*sin(Theta);

a = min(Theta(:)); b = max(Theta(:));

inttheta = integral(fint(n,m,R,Theta,Phi),a,b)

This works as I'd expect by assigning function handles appropriately, however when the integral line runs I recieve the following error:

Error using *

Arguments must be 2-D, or at least one argument must be scalar. Use TIMES (.*) for elementwise multiplication.

Running dbstop if error and playing with the workspace gives each section in fint evaluating appropriately (evaluated suggestion) i.e finnm, rYnm, sin(Theta) at this point will all run and give the NxNxN matrix I'd expect individually. I'm therefore stumped on what to do.

[Star Strider]

I’m lost. (I also don’t see where the ‘D’ matrix comes into all this, although I may have missed something.)

You need to use element-wise operations here:

        fint = @(n,m,R,Theta,Phi)...
            funnm(R,Theta,Phi).*rYnm(n,m,R,Theta,Phi).*sin(Theta);    

You may also need to use the ('ArrayValued',true) name-value pair, as I did in my example code.

Here, you need to specify the variable-of-integration:

        inttheta = integral(fint(n,m,R,Theta,Phi),a,b)

Since you defined ‘a’ and ‘b’ as the limits of ‘Theta’, this may be appropriate:

        inttheta = integral(@(Theta) fint(n,m,R,Theta,Phi), a, b)

noting that all the other variables and arguments must already be defined in your workspace.

If you are doing multiple integrations, you can use the related integral functions. However only integral accepts ('ArrayValued',true). If your multiple integral is array-valued, you need to use it with the integral function to create an appropriate multiple integration with array-valued evaluations integrating with respect to each variable.