Search code examples
matlabnumerical-integration

Matlab integration issue, issue with input function


I am trying to do a numerical integration on a function which is absolutely horrible to expand and work with analytically. The integral is over dpsi and dtheta. If I store the variables as sym i am told the integral input must be a double or single, if i store it as a double or single i am told I am adding two tensors of different dimensions. Help?

eta = input('Enter Dielectric Constant 1.5-4:  ');
psi = input('Enter Lattitude -pi/2 to +pi/2:  ');
theta = input('Enter Longitude -pi/2 to +pi/2:  ');
sdev = input('Enter STD DEV (roughness) maybe 0.1:  ');



dpsi = sym('dspi');
dtheta = sym('dtheta');

calpha = (cos(theta+dtheta)).*(cos(psi+dpsi));
rp01 = calpha-sqrt(eta-1+((calpha).^2));
rp02 = calpha+sqrt(eta-1+((calpha).^2));
rperp = (rp01./rp02).^2;
rp11 = ((eta.*calpha)-sqrt(eta-1+((calpha).^2)));
rp12 = ((eta.*calpha)+sqrt(eta-1+((calpha).^2)));
rpar = (rp11./rp12).^2;



fun = @(dtheta,dpsi) (rpar+rperp)
thetamax = (pi/2) - theta
psimax = (pi/2) - psi

q = integral2(fun,-pi/2,thetamax,-pi/2,psimax)

Solution

  • Simplest solution:

    Convert the symbolic expression to an actual numerical function using:

    fun = matlabFunction((rpar+rperp),'vars',{dtheta,dpsi});
    

    Better solution:

    Avoid the symbolic stuff in the first place and just define a function:

    function out = YOURFANCYFUNCTION(eta,psi,theta,sdev,dtheta,dpsi)
    calpha = (cos(theta+dtheta)).*(cos(psi+dpsi));
    rp01 = calpha-sqrt(eta-1+((calpha).^2));
    rp02 = calpha+sqrt(eta-1+((calpha).^2));
    rperp = (rp01./rp02).^2;
    rp11 = ((eta.*calpha)-sqrt(eta-1+((calpha).^2)));
    rp12 = ((eta.*calpha)+sqrt(eta-1+((calpha).^2)));
    rpar = (rp11./rp12).^2;
    out = (rpar+rperp);
    

    The function you pass to integral should then be:

    fun = @(dtheta,dpsi) YOURFANCYFUNCTION(eta,psi,theta,sdev,dtheta,dpsi)
    

    Which captures the current values of eta,psi,theta,sdev and makes dtheta and dpsi variables.

    By the way: The variable sdev is never used.