Error "Limits of integration must be double or single scalars." when solve equations with "fsolve" function

I want to solve the balance equations like following:

% momentum and force balance

f1 = int(pl*x,x,-L,0)+int(pr*x,x,0,L)+M/d;

f2 = int(pl,x,-L,0)*sin(alph+tilt)+int(pr,x,0,L)*sin(alph-tilt)-G/d-2*L*P_e*sin(alph)*cos(tilt);

f3 = int(pl,x,-L,0)*cos(alph+tilt)-int(pr,x,0,L)*cos(alph-tilt)+2*L*P_e*sin(alph)*sin(tilt);

initial_guess = [0.5, 0.1, 2];

equations = matlabFunction(f1,f2,f3, 'Vars', {x1, y1, Uratio});

result = fsolve(@(vars) equations(vars(1), vars(2), vars(3)), initial_guess);

Then error happend after running:

Error using integral

Limits of integration must be double or single scalars.

Error in symengine>@(x)x.*(Uratio.^3.*(Uratio.*(x.^4.*(-cos(pi./6.0-atan(x1./y1).*2.0).*(x1.^2./4.0+y1.^2./4.0)+cos(pi./6.0-atan(x1./y1).*2.0).*(x1.^2.*(3.0./4.0)+y1.^2.*(3.0./4.0))+x1.^2+y1.^2)+(x.^2.*(x1.^2+y1.^2).^2)./2.0+x.^6./6.0-x.^3.*cos(pi./1.2e+1-ata.........

Error in integralCalc/iterateScalarValued (line 314)

fx = FUN(t);

Error in integralCalc/vadapt (line 132)

[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)

[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 87)

Q = integralCalc(fun,a,b,opstruct);

Error in symengine>@(x1,y1,Uratio)deal(integral(@(x)x.*(Uratio.^3.*(Uratio.*(x.^4.*(-cos(pi./6.0-atan(x1./y1).*2.0).*(x1.^2./4.0+y1.^2./4.0)+cos(pi./6.0-atan(x1./y1).*2.0).*(x1.^2.*(3.0./4.0)+y1.^2.*(3.0./4.0))+x1.^2+y1.^2)+(x.^2.*(x1.^2+y1.^2).^2)./2.0+x.^6./6.0-x.^3.*cos(pi./1.2e+1-atan(x1./y1)).*(x1.^2+y1.^2).^(3.0./2.0).*(4.0./3.0)-x.^5.*cos(pi./1.2e+1-atan(x1./y1)).*sqrt(x1.^2+y1.^2).*(4.0./5.0))+(Uratio.*integral(@(x)(sin(pi./1.2e+1-atan(x1./y1)).^2.*(x1.^2+y1.^2)+(x-cos(pi./1.2e+1-atan(x1./y1)).*sqrt(x1.^2+y1.^2)).^2).^(3.0./2.0),0.0,x).*(cos(pi./6.0-atan(x1./y1).*2.0).*(x1.^2.*4.5265984e-5+y1.^2.*4.5265984e-5)-cos(pi.*(1.1e+1./6.0)-atan(x1./y1).*2.0).*(x1.^2.*4.5265984e-5+y1.^2.*4.5265984e-5)-cos(pi./6.0-atan(x1./y1).*2.0).*(x1.^2.*1.357.........

Error in fsolvetest_3para>@(vars)equations(vars(1),vars(2),vars(3)) (line 74)

result = fsolve(@(vars) equations(vars(1), vars(2), vars(3)), initial_guess);

Error in fsolve (line 264)

fuser = feval(funfcn{3},x,varargin{:});

Error in fsolvetest_3para (line 74)

result = fsolve(@(vars) equations(vars(1), vars(2), vars(3)), initial_guess);

Caused by:

Failure in initial objective function evaluation. FSOLVE cannot continue.

[Yuting]

Soory for the long error report

[Yuting]

the entire code is:

alph=pi*15/180;

L=0.116;

T_ice=18;

addweight=0.212;

T0=15;

dL=0.02;

%tilt=2.5*pi/180;

tilt=0;

heightofweight=12.5;

forcefromcable=0;

mu_L=0.001;

rho_S=920*0.95;

rho_L=1000;

Cp_s=2049.41;

h_m=334000+Cp_s*T_ice;

K_L=0.57;

P_e=102770;

% 101300+1000*9.8*0.15

g1=addweight*9.8;

g2=0.38*9.8;

heightofmasscenter=0.07;

heightofcable=0.13;

syms x1 y1 Uratio x

d = pi*L*sin(alph)/2;

M = g1*(dL+heightofweight*tan(tilt))*cos(tilt)+g2*heightofmasscenter*sin(tilt)-forcefromcable*heightofcable*sin(tilt);

G = g1+g2;

beta = atan(x1/y1);

xL = sqrt(x1^2+y1^2)*cos(pi-alph-beta);

DL = sqrt(x1^2+y1^2)*sin(pi-alph-beta);

xR = sqrt(x1^2+y1^2)*cos(beta-alph);

DR = sqrt(x1^2+y1^2)*sin(beta-alph);

LL = DL^2+(xL-x)^2;

RR = DR^2+(xR-x)^2;

C = Uratio*(int(x*RR^2,x,0,L)-int(x*LL^2,x,0,-L))/(int(LL^1.5,x,0,-L)-int(RR^1.5,x,0,L));

P0 = 12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*int(LL^2*x,x,0,-L)+C*int(LL^1.5,x,0,-L))/(rho_L*T0^3*K_L^3)+P_e;

pl = -12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*int(LL^2*x,x,0,x)+C*int(LL^1.5,x,0,x))/(rho_L*T0^3*K_L^3)+P0;

pr = -12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*int(RR^2*x,x,0,x)+C*int(RR^1.5,x,0,x))/(rho_L*T0^3*K_L^3)+P0;

f1 = int(pl*x,x,-L,0)+int(pr*x,x,0,L)+M/d;

f2 = int(pl,x,-L,0)*sin(alph+tilt)+int(pr,x,0,L)*sin(alph-tilt)-G/d-2*L*P_e*sin(alph)*cos(tilt);

f3 = int(pl,x,-L,0)*cos(alph+tilt)-int(pr,x,0,L)*cos(alph-tilt)+2*L*P_e*sin(alph)*sin(tilt);

disp('balance equations done')

initial_guess = [0.5, 0.1, 2];

equations = matlabFunction(f1,f2,f3, 'Vars', {x1, y1, Uratio});

result = fsolve(@(vars) equations(vars(1), vars(2), vars(3)), initial_guess);

Answer 1

[Harald]

Hi Yuting,

since you are solving the equations numerically anyway, my recommendation would be to perform all the computations numerically from the start, i.e. write a function that computes [f1, f2, f3] from (x1, y1, Uratio). I would only make an exception from this if the symbolic computations are really simple or can be simplified significantly.

This should not only help resolve the error but it will likely also improve performance.

Best wishes,

Harald

[Yuting]

Thanks for your reply!!

Do you mean I need to skip all the transfer variaties like xL, xR, DL, DR, LL, RR, C, P0, pl, pr... and calculate the f1, f2, f3 directly from the x1, y1, Uratio?

Why should I write a function about this?

Best,

Yuting

[Harald]

I would not skip them, but calculate them depending on the currently symbolic variables. You could highlight the lines between between "syms..." and "disp", then try the "Convert to Function" from the context menu to get started with such a function. You may need to switch some symbolic expressions to anonymous function handles and calls of int may need to be switched to integral.

At the moment, the formulas are just getting too long for me to even understand what exactly the problem is. This is the reason for me to suggest the workaround.

Best wishes,

Harald

[Yuting]

Hi @Harald,

Thanks for your suggestion! And we find a place which might cause the error when we rewrited the formulas.

In the following two lines, the limits in function "int" are 0 and x:

pl = -12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*int(LL^2*x,x,0,x)+C*int(LL^1.5,x,0,x))/(rho_L*T0^3*K_L^3)+P0;

pr = -12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*int(RR^2*x,x,0,x)+C*int(RR^1.5,x,0,x))/(rho_L*T0^3*K_L^3)+P0;

After replaced the int into integral, the upper limit x caused the error as it's not a scalar.

Do you have any possible solution?

Thanks again!

Yuting

[Harald]

Hi @Yuting,

I would avoid using the same variable as the integration variable and the limit as it has potential for confusion.

Since the integral will depend on the upper limit, I would create an anonymous function handle for this:

pl = @(x) -12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*integral(@(xi) LL(xi)^2*xi,0,x)+C*integral(@(xi) LL(xi)^1.5, 0,x))/(rho_L*T0^3*K_L^3)+P0;

This assumes that LL has been created as a function handle, e.g.

LL = @(x) DL^2+(xL-x)^2;

I hope this helps.

Best wishes,

Harald