Evaluating integrals with array limits

Jordan Schofield

2017 4 26

0 个回答

5 次查看（30 天）

1 个投票

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

Jordan Schofield 2017-4-27

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