Error: Limits of integration must be double or single scalars. While solving equations numerically with symbolic integral limits

Problem Description:

The equation system to be solved {f1=0, f2=0, f3=0} contains variables pl and pr. The expressions for pl and pr involve integrals with variables as their limits (indicated by bold parts in the code). After processing through matlabFunction to create anonymous functions, "int" is replaced with "integral". However, "integral" does not support variables as limits, resulting in an error:

Error using integral

Limits of integration must be double or single scalars.

Potentially useful reference:

https://www.zhihu.com/question/46915508

This is a similar problem I found, but I can't understand it...

Problem Code:

++++++++++++++++++++++++++

alph=pi*15/180;
L=0.116;
T_ice=18;
addweight=0.212;
T0=15;
dL=0.02;
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;
g1=addweight*9.8;
g2=0.38*9.8;
heightofmasscenter=0.07;
heightofcable=0.13;

syms x1 y1 x Uratio

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,0,x)+C*int(LL^1.5, 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,0,x)+C*int(RR^1.5, 0,x))/(rho_L*T0^3*K_L^3)+P0;

% Perform definite integrals on pl and pr to obtain the target equation system
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); 

% Proceed with solving
% Initial guess for variables
initial_guess = [0.5, 0.1, 2];
% Define solving function
equations = matlabFunction(f1,f2,f3, 'Vars', {x1, y1, Uratio});
% Use fsolve to solve the equation system
result = fsolve(@(vars) equations(vars(1), vars(2), vars(3)), initial_guess);

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-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))+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)).*5.715599983683794e+20-1.0)

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))+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)).*5.715599983683794e+20-1.0),0.0,2.9e+1./2.5e+2)+integral(@(x)-x.*(Uratio.^3.*(Uratio.*(x.^4.*(-cos(pi.*(1.1e+1./6.0)-atan(x1./y1).*2.0).*(x1.^2./4.0+y1.^2./4.0)+cos(pi.*(1.1e+1./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.1e+1./1.2e+1)-atan(x1./y1)).*(x1.^2+y1.^2).^(3.0./2.0).*(4.0./3.0)-x.^5.*cos(pi.*(1.1e+1./1.2e+1)-atan(x1./y1)).*sqrt(x1.^2+y1.^2).*(4.0./5.0))+integral(@(x)(sin(pi.*(1.1e+1./1.2e+1)-atan(x1./y1)).^2.*(x1.^2+y1.^2)+(x-cos(pi.*(1.1e+1./1.2e+1)-atan(x1./y1)).*sqrt(x1.^2+y1.^2)).^2).^(3.0./2.0),0.0,x)).*5.715599983683794e+20-1.0),-2.9e+1./2.5e+2,0.0)+8.810850563679611e-1,integral(@(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))+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)).*(-5.715599983683794e+20)+1.0,0.0,2.9e+1./2.5e+2).*2.588190451025207e-1+integral(@(x)Uratio.^3.*(Uratio.*(x.^4.*(-cos(pi.*(1.1e+1./6.0)-atan(x1./y1).*2.0).*(x1.^2./4.0+y1.^2./4.0)+cos(pi.*(1.1e+1./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.1e+1./1.2e+1)-atan(x1./y1)).*(x1.^2+y1.^2).^(3.0./2.0).*(4.0./3.0)-x.^5.*cos(pi.*(1.1e+1./1.2e+1)-atan(x1./y1)).*sqrt(x1.^2+y1.^2).*(4.0./5.0))+integral(@(x)(sin(pi.*(1.1e+1./1.2e+1)-atan(x1./y1)).^2.*(x1.^2+y1.^2)+(x-cos(pi.*(1.1e+1./1.2e+1)-atan(x1./y1)).*sqrt(x1.^2+y1.^2)).^2).^(3.0./2.0),0.0,x)).*(-5.715599983683794e+20)+1.0,-2.9e+1./2.5e+2,0.0).*2.588190451025207e-1-6.293948740487748e+3,integral(@(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))+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)).*(-5.715599983683794e+20)+1.0,0.0,2.9e+1./2.5e+2).*(-9.659258262890683e-1)+integral(@(x)Uratio.^3.*(Uratio.*(x.^4.*(-cos(pi.*(1.1e+1./6.0)-atan(x1./y1).*2.0).*(x1.^2./4.0+y1.^2./4.0)+cos(pi.*(1.1e+1./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.1e+1./1.2e+1)-atan(x1./y1)).*(x1.^2+y1.^2).^(3.0./2.0).*(4.0./3.0)-x.^5.*cos(pi.*(1.1e+1./1.2e+1)-atan(x1./y1)).*sqrt(x1.^2+y1.^2).*(4.0./5.0))+integral(@(x)(sin(pi.*(1.1e+1./1.2e+1)-atan(x1./y1)).^2.*(x1.^2+y1.^2)+(x-cos(pi.*(1.1e+1./1.2e+1)-atan(x1./y1)).*sqrt(x1.^2+y1.^2)).^2).^(3.0./2.0),0.0,x)).*(-5.715599983683794e+20)+1.0,-2.9e+1./2.5e+2,0.0).*9.659258262890683e-1)

Error in solution>@(vars)equations(vars(1),vars(2),vars(3)) (line 52)
result = fsolve(@(vars) equations(vars(1), vars(2), vars(3)), initial_guess);

Error in fsolve (line 267)
fuser = feval(funfcn{3},x,varargin{:});

Caused by:
Failure in initial objective function evaluation. FSOLVE cannot continue.

++++++++++++++++++++++++++

[Dyuman Joshi]

Please attach the parameter definition section as well.

[Yuting]

@Dyuman Joshi Hi, the parameters were attached~

[Torsten]

C and P0 are missing.

The limit of the integrals

pl = ...int(LL^2*x,0,x)+C*int(LL^1.5, 0,x)...;
pr = ...int(RR^2*x,0,x)+C*int(RR^1.5, 0,x)...;

must be a fixed double, not a symbolic variable. Or is the limit of integration also a solution variable which should appear as "Vars" in the call to "matlabFunction" ?

[Yuting]

@Torsten Hi, the C and P0 were added.

And I have tried:

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

Sadly, another error accoured:

Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters.

[Torsten]

I don't know where this error message comes from, but setting

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

would mean that you want to solve for 4 variables with 3 equations. This will not be possible with "fsolve".

[Walter Roberson]

Solving 3 equations in 4 variables will not typically fail in fsolve. If it were the other way around, more variables than equations, that would typically fail

[Torsten]

Ok. "fsolve" might be able to cope with the problem, but usually there is something wrong with these tasks :-)

[Torsten]

@Yuting

Maybe you have double integrals that you have to compute where the inner integrals give functions of x, and x is integrated out by the outer integration ?

[Yuting]

@Torsten Yes, the pl and pr give functions of x as the inner integrals, and the x is integrated out by f1-f3 as outer integrals.

[Torsten]

In this case, set up your complete problem numerically with no symbolic variables and "integral2".

initial_guess = [0.5, 0.1, 2];
% Use fsolve to solve the equation system
result = fsolve(@(vars) fun(vars(1), vars(2), vars(3)), initial_guess);

function res = fun(x1,y1,Uratio)
  alph=pi*15/180;
  L=0.116;
  T_ice=18;
  addweight=0.212;
  T0=15;
  dL=0.02;
  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;
  g1=addweight*9.8;
  g2=0.38*9.8;
  heightofmasscenter=0.07;
  heightofcable=0.13;
  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 = @(x)DL^2+(xL-x).^2;
  RR = @(x)DR^2+(xR-x).^2;
  ...
  res = [f1;f2;f3]
end

[Yuting]

@Torsten Thanks！！！

I rewrote the codes like this, but there should be some mistakes... Could you help me to have a roughly check?

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;

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 = @(x1,y1) atan(x1/y1);
xL = @(x1,y1) sqrt(x1^2+y1^2)*cos(pi-alph-beta);
DL = @(x1,y1) sqrt(x1^2+y1^2)*sin(pi-alph-beta);
xR = @(x1,y1) sqrt(x1^2+y1^2)*cos(beta-alph);
DR = @(x1,y1) sqrt(x1^2+y1^2)*sin(beta-alph);

LL = @(x) DL^2+(xL-x)^2;
RR = @(x) DR^2+(xR-x)^2;
C = @(Uratio,x) Uratio*(integral(@(x)x*RR(x)^2,0,L)-integral(@(x)x*LL(x)^2,0,-L))/(integral(LL(x)^1.5,0,-L)-integral(RR(x)^1.5,0,L));
P0 = @(Uratio,x) 12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*integral(@(x)LL(x)^2*x,0,-L)+C*integral(LL(x)^1.5,0,-L))/(rho_L*T0^3*K_L^3)+P_e;
pl = @(Uratio,x) -12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*LL(x)^2*x+C*LL(x)^1.5)/(rho_L*T0^3*K_L^3)+P0;
pr = @(Uratio,x) -12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*RR(x)^2*x+C*RR(x)^1.5)/(rho_L*T0^3*K_L^3)+P0;

f1 = @(x1,y1,Uratio,x) integral2(@(xi,x)pl(Uratio,xi)*x,0,x,-L,0)+integral2(@(xi,x)pr(Uratio,xi)*x,0,x,0,L)+M/d; 
f2 = @(x1,y1,Uratio,x) integral2(@(xi,x)pl(Uratio,xi),0,x,-L,0)*sin(alph+tilt)+integral2(@(xi,x)pr(Uratio,xi),0,x,0,L)*sin(alph-tilt)-G/d-2*L*P_e*sin(alph)*cos(tilt); 
f3 = @(x1,y1,Uratio,x) integral2(@(xi,x)pl(Uratio,xi),0,x,-L,0)*cos(alph+tilt)-integral2(@(xi,x)pr(Uratio,xi),0,x,0,L)*cos(alph-tilt)+2*L*P_e*sin(alph)*sin(tilt); 
disp('balance equations done')

fun = {@(vars) f1(vars(1), vars(2), vars(3)); @(vars) f2(vars(1), vars(2), vars(3)); @(vars) f3(vars(1), vars(2), vars(3))};

initial_guess = [0.5, 0.1, 2];
result = fsolve(@(vars) cellfun(@(f) f(vars), fun), initial_guess);

[Torsten]

I still cannot understand your integrals.

Assume for a moment

LL = DL^2+(xL-x)^2;
pl = int(LL^2*x,0,x)
f1 = int(pl*x,x,-L,0)

How would f1 look like in a mathematical notation as a double integral ?

[Walter Roberson]

pl and pr have int() calls that are missing the variable of integration. The variable that will be assumed is not necessarily going to be the one you want.

[Yuting]

@Torsten It would be something like this:

or using integral2() with xi as another variable needs to be integrated:

[Yuting]

@Walter Roberson Thank you for your help!! I added every variables to every anonymous function as you suggested, but the error "Not enough input arguments." still exists.

beta = @(x1,y1) atan(x1/y1);

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

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

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

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

% 积分常数项

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

RR = @(x,x1,y1) DR^2+(xR-x)^2;

C = @(Uratio,x,x1,y1) Uratio*(integral(@(x)x*RR(x,x1,y1)^2,0,L)-integral(@(x)x*LL(x,x1,y1)^2,0,-L))/(integral(LL(x,x1,y1)^1.5,0,-L)-integral(RR(x,x1,y1)^1.5,0,L));

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

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

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

% momentum balance

f1 = @(x1,y1,Uratio,x) integral2(@(xi,x)pl(Uratio,xi,x1,y1)*x,0,x,-L,0)+integral2(@(xi,x)pr(Uratio,xi,x1,y1)*x,0,x,0,L)+M/d;

% Force balance分解为竖+横，左+右

f2 = @(x1,y1,Uratio,x) integral2(@(xi,x)pl(Uratio,xi,x1,y1),0,x,-L,0)*sin(alph+tilt)+integral2(@(xi,x)pr(Uratio,xi,x1,y1),0,x,0,L)*sin(alph-tilt)-G/d-2*L*P_e*sin(alph)*cos(tilt);

f3 = @(x1,y1,Uratio,x) integral2(@(xi,x)pl(Uratio,xi,x1,y1),0,x,-L,0)*cos(alph+tilt)-integral2(@(xi,x)pr(Uratio,xi,x1,y1),0,x,0,L)*cos(alph-tilt)+2*L*P_e*sin(alph)*sin(tilt);

disp('balance equations done')

fun = {@(vars) f1(vars(1), vars(2), vars(3)); @(vars) f2(vars(1), vars(2), vars(3)); @(vars) f3(vars(1), vars(2), vars(3))};

Answer 1

[Torsten]

I guess this is what you mean.

initial_guess = [0.5, 0.1, 2];
% Use fsolve to solve the equation system
result = fsolve(@(vars) fun(vars(1), vars(2), vars(3)), initial_guess, optimset('MaxFunEvals',10000,'MaxIter',10000));

function res = fun(x1,y1,Uratio)
  alph=pi*15/180;
  L=0.116;
  T_ice=18;
  addweight=0.212;
  T0=15;
  dL=0.02;
  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;
  g1=addweight*9.8;
  g2=0.38*9.8;
  heightofmasscenter=0.07;
  heightofcable=0.13;
  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 = @(x)DL^2+(xL-x).^2;
  RR = @(x)DR^2+(xR-x).^2;
  
  C = Uratio*(integral(@(x)x.*RR(x).^2,0,L)-integral(@(x)x.*LL(x).^2,0,-L))/(integral(@(x)LL(x).^1.5,0,-L)-integral(@(x)RR(x).^1.5,0,L));
  P0 = 12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*integral(@(x)LL(x).^2.*x,0,-L)+C*integral(@(x)LL(x).^1.5,0,-L))/(rho_L*T0^3*K_L^3)+P_e;
  pl = @(y)(-12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*integral(@(x)LL(x).^2.*x,0,y)+C*integral(@(x)LL(x).^1.5, 0,y))/(rho_L*T0^3*K_L^3)+P0);
  pr = @(y)(-12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*integral(@(x)RR(x).^2.*x,0,y)+C*integral(@(x)RR(x).^1.5, 0,y))/(rho_L*T0^3*K_L^3)+P0);
  % Perform definite integrals on pl and pr to obtain the target equation system
  f1 = integral(@(y)pl(y).*y,-L,0,'ArrayValued',1)+integral(@(y)pr(y).*y,0,L,'ArrayValued',1)+M/d; 
  f2 = integral(@(y)pl(y),-L,0,'ArrayValued',1)*sin(alph+tilt)+integral(@(y)pr(y),0,L,'ArrayValued',1)*sin(alph-tilt)-G/d-2*L*P_e*sin(alph)*cos(tilt); 
  f3 = integral(@(y)pl(y),-L,0,'ArrayValued',1)*cos(alph+tilt)-integral(@(y)pr(y),0,L,'ArrayValued',1)*cos(alph-tilt)+2*L*P_e*sin(alph)*sin(tilt); 
  res = [f1;f2;f3];
end

[Yuting]

Thanks!! The code is running without any error report now! Now my mission is to find the right initial guess values

Thanks again!

[Torsten]

... and search for possible errors in your equations.

E.g. the fact that some lower limits in your integrals are bigger than the upper limits looks suspicious.

[Yuting]

Okay! I will definitely do it

Thanks again!

Best wishes,

Yuting

[Yuting]

Hi @Torsten!

Sorry to bother you again!

I found it's kind of hard to search the right initial guess value manually is difficult... Should I write a loop to sweep the target area? Or is there any more convinient method?

Best,

Yuting

[Torsten]

I don't know if "fsolve" can be combined with "MultiStart".

If not, reformulate your problem as

result = fmincon(@(vars) fun(vars(1), vars(2), vars(3)), initial_guess, [],[],[],[],[],[],[],optimset('MaxFunEvals',10000,'MaxIter',10000));

instead of

result = fsolve(@(vars) fun(vars(1), vars(2), vars(3)), initial_guess, optimset('MaxFunEvals',10000,'MaxIter',10000));

and

res = f1^2+f2^2+f3^3   

instead of

res = [f1;f2;f3];

"fmincon" accepts "MultiStart" - if possible with lower and upper bounds on the solution variables.

Instead of "fmincon", you can also try "lsqnonlin" as solver in the new formulation.

[Yuting]

Thanks!!!

[Torsten]

Sorry - lsqnonlin must be used in the old formulation:

result = lsqnonlin(@(vars) fun(vars(1), vars(2), vars(3)), initial_guess, [],[], optimset('MaxFunEvals',10000,'MaxIter',10000));

and

res = [f1;f2;f3];

Answer 2

[Walter Roberson]

integral() can never accept symbolic integration limits.

If you are trying to build up a symbolic algorithm that is later to be used with matlabFunction(), then code the integral() as either int() or vpaintegral()

Note: integral() does not accept expressions such as pl*x . The first parameter to integral() must be a function handle.

[Yuting]

Thanks! I have fixed the wrong usage of "integral" in the code~

And if the symbolic limits are never accepted by integral(), is there any other possibilities to solve the equations numerically?

I can't expand this part int(LL^2*x,0,x)+C*int(LL^1.5, 0,x) into polynomial because of the ^1.5.

Best,

Yuting

[Walter Roberson]

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

Notice that you have x as an upper bound there.

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

and that int() with x as an upper bound is nested within a second int over x.

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

when matlabFunction() converts those int() into calls to integral(), it is not especially smart about it. It creates something with roughly the form

integral(@(x) expression+integral(@(x)second_expression, 0, x), 0, 0.116)

so one integral() is calling another, and the second one is using the variable of the first as the upper bound.

However, by default, integral() passes in a vector of values -- so x is a vector at the time it is received by the inner integral() call that tries to use it as the upper bound. That's a problem since the bounds can only be scalar.

What can you do?

Well, if you were tell matlabFunction to write to a 'file' and were to edit the code slightly, you could edit the outer integral() call to specify 'ArrayValued', true as its options. That would cause it to pass scalar x to the anonymous function, and so the upper bound of the inner integral() would end up being a scalar.

But to fix this automatically without doing that kind of editting... is a bit of a nuisance.

[Walter Roberson]

Note by the way you have

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

but that should be

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

Also I suggest you consider using

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

as that would allow you to

result = fsolve(equations, initial_guess);

[Yuting]

Thanks a lot for your reply!!

The update of "equations = matlabFunction([f1,f2,f3], 'Vars', {[x1, y1, Uratio]})" works!

But I'm still confused with the integral() stuff...

I wrote something like this:

matlabFunction(expression + integral(@(x) second_expression, 0, x), 'File', 'outer_integral_function', 'ArrayValued', true);

I hope this is correct and I don't know what to do next at all... Sorry for my poor matlab skill...