Vapsolve fails using limits

Question

Jesus Mezquita 2018-6-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/405311-vapsolve-fails-using-limits

评论： Jesus Mezquita 2018-6-14

I am working in a script to solve two equations with two unknowns. I have an error when I establish the limits for the variable Qv and Pw. The outputs values for Qv should be only positivies and the output values for Pw should be between 0 and 1 in each iteration. If I don't use limits to constraint the result I finsih with negatives and complex numbers which is a solutions but doesn't work for my porpuse. I want to continue running the script avoiding the error or maybe use another function to get a result between these limits.

%-----Script1-----%
clear; clc; close all;
%---Import data
Data = importdata('ERT_It5.prn',' ',1)
%---Variables
Data = Data.data;
L = length(Data);   % Length of data variables
Rt = Data(:,3)'; 
a = Data(:,7)';  
m = Data(:,6)';  
n = Data(:,5)';  
Rw = Data(:,8)'; 
Sw = Data(:,4)'; 
    for x = 1:L
        P(x) = (((a(x)*Rw(x))/(Rt(x)*(Sw(x).^n(x)))).^(1/m(x))); 
    end
min(P),max(P) % Maximum and minimun values
figure(1)
histogram(P)
CEC = Data(:,9)';
d = Data(:,10)';
syms Qv Pw
assume(Qv,'real')
assume(Pw,'real')
    for y = 1:L
        B(y) = 0.046*(1-0.6*exp(-(1/Rw(y))/0.013));
        Eq1 = Qv == d(y)*(((1-Pw)/Pw))*CEC(y);  % Equation 1  
        Eq2 = Pw == ((a(y)*Rw(y)*Sw(y))/(1+B(y)*Qv*Rw(y)*Rt(y)*(Sw(y).^n(y)))).^(1/m(y));   % Equation 2
        [RQv(y),RPw(y)]= vpasolve(Eq1,Eq2,Qv,Pw,[0 Inf;0 1])
    end
    OUT = double(real(RPw'));
    Excel = arrayfun(@(x) sprintf('%10.6f',x),OUT,'un',0);
    xlswrite('Outputfile.xlsx',Excel)
%-----End-----%

ERROR

In an assignment  A(:) = B, the number of elements in A and B must be the
same.
Error in sym/privsubsasgn (line 1067)
                L_tilde2 =
                builtin('subsasgn',L_tilde,struct('type','()','subs',{varargin}),R_tilde);
                Error in sym/subsasgn (line 904)
            C = privsubsasgn(L,R,inds{:});
Error in script (line 34)
        [RQv(y),RPw(y)]= vpasolve(Eq1,Eq2,Qv,Pw,[0 Inf;0 1])