Evaluating integrals with array limits

Dear Matlab users,

I'm hope someone can help me with solving these integrals more efficiently. I have written a code to iteratively solve the function f(Sw) shown below, however, the values of Sw, Sw* and Swi (the limits of the integral) must all be single values instead of arrays for it to work. So, to create curves of f(Sw) vs. Sw, I need to run the function again and again while manually changing the value of Sw.

Troubleshooting so far: 1. Using the integral Matlab function with array limits, which returned the error "First input argument must be a function handle.". 2. Next, I made the functions within the integrals "Anonymous" and use the integral Matlab function, this returns the error "A and B must be floating-point scalars." 3. After, I tried using the int function which returned the error "Invalid argument."

So for now I am using the int function (as shown below) with singular values of Sw, Sw*, and Swi. The resulting curves look reasonable, but take a long time to generate.

clear; clc;
%%%%***Code for Relative Imbibition Rate Function*** %%%%
%%%Symbolic variables
syms Swv; syms Siwv; syms Swiv; syms Swmv; syms a;

%%%Water Saturations
Sw =    0.001;                  % Water saturation
Swi =   0.0;                    % Initial water saturation
Siw =   0.3180;                 % Minimum water saturation
Swm =   0.9850;                 % Maximum water saturation

%%%Relative Permeability
Krw = abs(((1-(1-(((Sw-Siw)./(Swm-Siw)).^1.49)).^0.671).^2).*(((Sw-Siw)./(Swm-Siw)).^0.5));

%%%Capillary Pressure (MPa)
Pc = 0.0872.*((((Swv-Swiv)./(Swmv-Swiv)).^-1.49)-1).^0.329;
dPc = diff(Pc,Swv);
dPc = subs(dPc,Swv,Sw);
dPc = subs(dPc,Swmv,Swm);
dPc = subs(dPc,Swiv,Swi);
dPc = [zeros(size(dPc,1),1) dPc];

%%%Integration variables
fSw =   ((Sw-Swi)/(Swm-Swi));      % Initial Guess of fSw
fSw_old = 0;
countf = 0;

%%Iterative Integration to find f(Sw)
while max(abs(fSw - fSw_old)) > 1e-8 

fSw_old = fSw;

Sw_1 = @(a) ((Sw-a).*Krw.*Pc)./fSw;
Sw_2 = @(a) ((Swi-a).*Krw.*Pc)./fSw;

fSw1 =  (int(Sw_1,a,Sw,Swm));
fSw2 =  (int(Sw_2,a,Swi,Swm));
fSw  =  1 - (fSw1./fSw2);

       countf = countf + 1;

      %Generate error if there are more than 100 iterations
       if countf > 100   
       error('check code'); 
       end 

end

fSw(isnan(fSw)==1) = 0 ;
fSw(isinf(fSw)==1) = 0 ;
disp(fSw);

Please kindly let me know of any troubleshooting ideas you may have if you're interested in the problem. Specifically, how do I solve these integrals with array limits?

Many thanks and kind regards, J

[Torsten]

I don't understand why you have to iterate for f(S_w).

Since you can take f(S_w) out of the two integrals as a constant value, it cancels out on the right-hand side of your equation.

Best wishes

Torsten.

[Jordan Schofield]

Hi Torsten,

Thanks for your comment, indeed you're right, I made an error in the integral denominators. I've corrected it above, so the integrals are integrated with respect to α, and the denominators are f(α). Given that, they can't be taken out of the integrals and cancelled.

Kind regards, Jordan

[Torsten]

Then I also don't see a reason to iterate.

And the f's in f(S_w) and f(alpha) are different - thus f(S_w) and g(alpha) ?

Best wishes

Torsten.