Keep getting Empty sym: 0-by-1 when attempting to solve bounds of an integral

Hi there,

I'm attempting to solve the equation for h when

W=r*b*int(Taux*sin(Pheta)+Sig*cos(Pheta),Pheta,[PhetaR,PhetaF])

Where:

PhetaR=acos(1-k*h/r)
PhetaF=acos(1-h/r)
Taux=((c1+(Kc1/b+Kphi1)*h^n1*tand(phi1))*(1-exp(-r*(PhetaF-Pheta-(1-s)*(sin(PhetaF)-sind(Pheta)))/kx)))
Sig=(Kc1/b+Kphi1)*h^n1

Below is the script I am using:

clc; clear all; close all;

r=0.1; 
b=0.15; 
Kc1=0.99*1000;
Kphi1=1528.4*1000;
n1=1.1;   
c1=1.04*1000;    
phi1=0.489;
k=0.8;
s=0.5;
kx=0.03'
W=100 

syms h real positive
syms Pheta real positive

PhetaR=acos(1-k*h/r)
PhetaF=acos(1-h/r)

fun1=((c1+(Kc1/b+Kphi1)*h^n1*tan(phi1))*(1-exp(-r*(PhetaF-Pheta-(1-s)*(sin(PhetaF)-sin(Pheta)))/kx)))*sin(Pheta)+((Kc1/b+Kphi1)*h^n1)*cos(Pheta)

eqn1=w==r*b*int(fun1,Pheta,[PhetaR,PhetaF])
 
vpasolve(eqn1,h,.1)

[Walter Roberson]

try changing the int() to vpaintegral

[Louis McDermott]

I've tried vpaintegral and it spits out a complex number, and I know the expected value should be real.

[Walter Roberson]

keep vpaintegral and

vpasolve(eqn1,h,[-inf, inf] )

You do not get your choice of starting point within the interval but you do get to confine to real.

[Louis McDermott]

Heres what i've done:

clc; clear all; close all;

r=0.1; 
b=0.15; 
Kc1=0.99*1000;
Kphi1=1528.4*1000;
n1=1.1;   
c1=1.04*1000;    
phi1=28;
k=0.8;
s=0.5;
kx=0.03'
W=200

Sig=(Kc1/b+Kphi1)
    
syms h real positive
syms Pheta real positive
    
fun1=(((c1+Sig*h^n1*tand(phi1))*(1-exp(-(r*(acos(1-h/r)-Pheta-(1-s)*(sin(acos(1-h/r))-sin(Pheta))))/kx)))*sin(Pheta)+(Sig*h^n1*cos(Pheta)))

eqn1=W==r*b*vpaintegral(fun1,Pheta,[acos(1-k*h/r),acos(1-h/r)])

vpasolve(eqn1,h,[-inf,inf])

However, after pressing run, it takes too long to solve (I left it running for 20 minutes and it couldnt solve for h).

[Walter Roberson]

Could you explain why you take the complex congugate transpose of the real scalar 0.03 for kx ?

Answer 1

[Walter Roberson]

Provided that you entered your coefficients and equations correctly (I did not check), there does not appear to be any real-valued solution.

Q = @(v) sym(v);
r = Q(0.1); 
b = Q(0.15); 
Kc1 = Q(0.99)*1000;
Kphi1 = Q(1528.4)*1000;
n1 = Q(1.1);   
c1 = Q(1.04)*1000;    
phi1 = Q(28);
k = Q(0.8);
s = Q(0.5);
kx = Q(0.03)' ;
W = Q(200);

Sig = (Kc1/b+Kphi1);
    
syms h real positive
syms Pheta real positive
    
fun1 = (((c1+Sig*h^n1*tand(phi1))*(1-exp(-(r*(acos(1-h/r)-Pheta-(1-s)*(sin(acos(1-h/r))-sin(Pheta))))/kx)))*sin(Pheta)+(Sig*h^n1*cos(Pheta)));
disp(fun1)

eqn1 = W==r*b*vpaintegral(fun1,Pheta,[acos(1-k*h/r),acos(1-h/r)]);
eqn1 = rhs(eqn1) - lhs(eqn1);

hvals = linspace(-0.01,0.2);
EQN1 = arrayfun(@(H) double(subs(eqn1, h, H)), hvals);
plot(hvals, real(EQN1), 'b', hvals, imag(EQN1), 'r');

Blue is real component, red is imaginary component. If you use data cursor you can see that to the left edge of the red, the imaginary starts slightly below 0, and then it runs at 0 for a while, eventually again becoming non-zero a little below 0.2. Any real-valued solution would have to have an imaginary component of 0, so only in the range plotted above. If you plot wider, the values quickly go away from zero.

But then when we plot the real component, the line in blue, we can see that it is decidedly not crossing zero anywhere in that range.

We should be able to narrow down the range for the function being real-valued. The sqrt(1-(10*h-1)^2) requires that 1-(10*h-1)^2 is non-negative, so 1-(10*h-1)^2 = 0 is a boundary. SO 1 = (10*h-1)^2 and sqrt() of both sides gives 1 = (10*h-1) which is 2 = 10*h -> h = 2/10 on one side. On the other side, -1 = 10*h-1 -> 0 = 10*h -> h = 0. So the function is only "naturally" real-valued between 0 and 0.2 . However, if 1-10*h is outside the range +/-1 then acos(1-10*h) could hypothetically be complex-valued in just the right amount to cancel the complex component of the sqrt(), and negative h would lead to complex h^1.1 so hypothetically there could be individually points (not runs) where the function had imaginary component 0; that would have to be chased down more rigourously. It would have to be outside any range that I tested. And it seems pretty unlikely that in such a case, the real components would also just happen to evaluate to 0. And for h beyond roughly +/-10^4, the expression overflows symbolic computation range.