Undefined function 'minus' for input arguments of type 'function_handle'. Error in flow1d2phase (line 46) Sw_shock = fzero(@(Sw​)(dfw(Sw)-​(fw(Sw0)-f​w(Sw))./(S​w0-Sw)), 0.5);

close;
clear;
clc;
 
%Reservoir properties
%ENTER reservoir properties
k = 300;
poro = 0.2;
A = 100;
miuo = 1;
miuw = 1;
L = 3000;
q = 1;
Bo = 1;
Bw = 1;
 
%ENTER IC and BCs
Pr = 29.007547546;
Pinit = 3000;
 
%Saturation
Sw0 = 0.2;
Swc = 0.2;
Sor = 0.15;
 
%Grid and time
%ENTER number of grid and time
t = [100, 200, 300, 400, 500, 600, 700];
 
nt = length(t);
 
%Functions
 
Se = @(Sw)((Sw-Swc)/(1-Swc-Sor));
dSe = 1/(1-Swc-Sor);
krw = @(Sw)0.5*(Se).^(6)+0.03*Se;
dkrw = @(Sw)3*(Se).^(5)+0.03;
kro = @(Sw)0.7*(1-Se).^(5)+0.05*(1-Se);
dkrw = @(Sw)-3.5*(1-Se).^(4)-0.95;
fw = @(Sw)((krw(Sw)/miuw)./(krw(Sw)/miuw+kro(Sw)/miuo));
dfw = @(Sw)((dkrw(Sw)/miuw.*kro(Sw)/miuo-krw(Sw)...
/miuw.*dkro(sw)/miuo)./(krw(sw)/miuw+kro(sw)/miuo).^2);
 
Sw0 = Swc;
 
Sw_shock = fzero(@(Sw)(dfw(Sw)-(fw(Sw0)-fw(Sw))./(Sw0-Sw)), 0.5);
 
Sw_inj = 1-sor;
 
Sw_rf = linspace(Sw_inj, Sw_shock, 100);
Sw_all = linspace(Swc, 1-Sor, 100);
 
xSw_all = (q*t/A*poro)*dfw(Sw_all);
 
xSw_rf= (q*t/A*poro)*dfw(Sw_rf);
xt_shock = xSw_rf(end);
 
figure(1);
subplot(3,1,2); plot(Sw_all, fw(Sw_all)); axis([0 1 0 1]);
xlabel('Sw'); ylabel('fw');
subplot(3,1,3);plot([xSw_rf xt_shock max(xSw_all)], ...
[Sw_rf Sw0 Sw0], 'r', xSw_all, Sw_all, '--b');
xlabel('x/t'); ylabel('Sw'); legend('BL solution');
 

Answer 1

[Walter Roberson]

Se = @(Sw)((Sw-Swc)/(1-Swc-Sor));

krw = @(Sw)0.5*(Se).^(6)+0.03*Se;

In krw, Se is not a named parameter to the anonymous function, so MATLAB is going to look through the workspace to see if it can find a variable named Se. It succeeds, finding the function handle you defined a few lines earlier. It records the content of the variable Se inside the anonymous function being built for krw

At the time that krw is invoked, it will be passed one parameter positionally. Anything directly in the body that is a match for that positional parameter's name, Sw will have the value of the positional parameter substituted when invoking the function. But there is nothing in the function body that is named Sw so no substitution will be done. The expression for the function body will then start executing. 0.5 is okay. Then the (Se) part. Se will not be dynamically looked up to find what current value Se might have in the workspace: instead, the value of Se that was recorded at the time the function handle was built will be recalled. So the handle to the anonymous function will be pulled out of the data structure representing krw .

Now you have 0.5 times an anonymous function handle. But you cannot do mathematics on function handles, so you get an error message.

If you wanted Se to be invoked on the passed in Sw then you needed to define krw differently:

krw = @(Sw)0.5*(Se(Sw)).^(6)+0.03*Se(Sw);

Anonymous function handle Se would be pulled out of the krw data structure, and would be invoked with the content of the positional parameter that was passed in to krw, which would cause a numeric result to be returned, and it is fine to raise the numeric result to a power and multiply the result by a numeric constant.

[Antonio Gomes]

Thanks for ur prompt reply!

Answer 2

[Steven Lord]

Se = @(Sw)((Sw-Swc)/(1-Swc-Sor));
dSe = 1/(1-Swc-Sor);
krw = @(Sw)0.5*(Se).^(6)+0.03*Se;

There are a couple places where you use Se (the anonymous function) on its own. You can't perform arithmetic like krw is trying to do (0.03*Se) on a function handle. You can evaluate that function handle and perform arithmetic on the result. One example is:

krw = @(Sw) 0.5*Se(Sw).^6 + 0.03*Se(Sw);

I'll leave it to you to modify the rest of the places your code tries to perform arithmetic on function handles.

[Antonio Gomes]

Thanks for ur prompt reply!