MATLAB: Discontinuities when computing integration of error functions using integral function

numerical integration

I am trying to integrate a function over a region in different time intervals. The integration looks something like this.

fun_uz = @(u)1./sqrt(u).*exp(-Z.^2./(2.*u));
fun_Y = @(u)(erf((Y+B)./sqrt(2.*u))-erf((Y-B)./sqrt(2.*u)));
fun_Z = @(u)(erf((X+L+u)./sqrt(2.*u))-erf((X-L+u)./sqrt(2.*u)));
fun = @(u)inc.*fun_uz(u).*fun_Y(u).*fun_Z(u);
fint = integral(fun,0,upperl);

The variable 'upperl' is the upper limit of the integral function. I have to perform this integration over different X,Y, and Z regions and different 'upperl' values. I am getting profiles which are discontinuous for different 'upperl' values. I have shown here profiles at few different 'upperl' values.

I am not able to understand why the discontinuity are occuring, any help is greatly appreciated. Thanks.

Best Answer

The function you are trying to integrate looks like following

Put these lines inside for loops

ezplot(fun,[0 upperl])
pause(0.01)

When time > 0.1 upperl is big. When you call integral you don't know how many points it takes. Maybe it misses something

Related Solutions

MATLAB: Double integral where one limit can only be solved numerically

The following is approximately what I think you need to do. It may need some adjusting here and there. I am unable to test it directly on my own primitive version of matlab, but at least it will give you an idea of what I was trying to say in the outline. Note that I have assumed that the value zero would be all right as a lower limit in the "estimate" interval passed to 'fzero'. It should be okay if p-q>=1. Otherwise setting it properly would be a more difficult task. The upper end of the interval, 'up', is doubled until the inequality reverses. I also assume r>=0 to ensure that p-q-lambda-sqrt(r*lambda+exp(lambda-x) decreases as lambda increases. I am assuming the parameters p,q, r, and V can be passed to the subfunction 'samymax' in the manner I've shown here. Let me know if you encounter unsurmountable difficulties with this. Note that I haven't placed any options as to error tolerances in either 'fzero' or 'integral2'. You may want to specify these to guarantee good accuracy.

I do have one question. Do you anticipate encountering any singularities in the integrand from the division by f(x,y) for the specified ranges of x and y? Matlab's 'integral2' function is supposed to be able to handle singularities in the integrand but only if they are of an integrable nature- that is, only if the integral isn't divergent.

Here is my suggested code:

 % The parameters
 p = ... % We assume p-q >= 1
 q = ...
 r = ...
 V = ...
 % The integrand as an anonymous function of x and y
 samint = @(x,y)A*exp(y-x)./f(x,y); % <-- Fix this up so it's correct
 % The inner integral upper limit as a function of x
 function lambda = samymax(X)
   [m,n] = size(X);
   lambda = zeros(m,n);
   for ix1 = 1:m
     for ix2 = 1:n
       x = X(ix1;ix2);
       up = 1; % We assume r >= 0
       while p-q-up-sqrt(r*up+exp(up-x) >= 0
         up = 2*up;  % Double 'up' each time until inequality is "<"
       end
       pqrx = @(lam)p-q-lam-sqrt(r*lam+exp(lam-x);
       % Assume p-q>=1, so that lower interval value of zero is valid 
       lambda(ix1,ix2) = fzero(pqrx,[0,up]);
     end
   end
 end
 % The double integral itself
 I = integral2(samint,0,V,0,@samymax);

MATLAB: Integral: Is it better to add integrals with different bounds or use one integral with a sum of characteristic functions as the integrand

Should you use Heaviside step functions to create an integral that is a piecewise discontinuous function, NO. NO. NO!!!!!! There, I said it three times, so it must be true.

Never intentionally create a function with discontinuities and then try to integrate it. Even if integral succeeds, the integration will be far more inefficient than it would be by simply splitting the problem up into segments. Use of a step function to combine pieces together is simply a lazy trick, hoping that integral will be smart enough to survive. Never be tricky, and hope that a computer algorithm will survive your trickery. Too often, it won't be.

Note that integral has the capability to provide waypoints. If you know there is a problem in some spot, then you need a waypoint.

Problems might be a discontinuity, or a singularity of some sort, or something that would cause the integration to fail. For example, the presence of a near delta function.

fun = @(x) exp(-50*(x-31.25).^2);
integral(fun,-1000,1000)
ans =
      0
integral(fun,-1000,1000,'waypoints',[30 31.25 32])
ans =
      0.25066

Yes, a waypoint will be a great improvement. For example, if we tried to integrate a simple step function.

step = @(x) double (x>pi);
timeit(@() integral(step,0,12))
ans =
    0.0020337
timeit(@() integral(step,0,12,'waypoints',pi))
ans =
   0.00052395
format long g
integral(step,0,12)
ans =
          8.85840830988541
integral(step,0,12,'waypoints',pi)
ans =
          8.85840734641021
12 - pi
ans =
          8.85840734641021

As you can see, the form with a waypoint is both significantly faster, as well as wildly more accurate.

Best Answer

Related Solutions

MATLAB: Double integral where one limit can only be solved numerically

MATLAB: Integral: Is it better to add integrals with different bounds or use one integral with a sum of characteristic functions as the integrand

Related Question