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

2023 8 6
2 个回答
回答已采纳
更新时间：2023 8 11
5 次查看（30 天）
Ran in:

0 个投票

Ran in:
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.
++++++++++++++++++++++++++

15 个评论
显示 13更早的评论 隐藏 13更早的评论

Torsten
Torsten 2023-8-6
编辑：Torsten 2023-8-6
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
Yuting 2023-8-6
@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
Torsten 2023-8-6
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
Walter Roberson 2023-8-6
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
Torsten 2023-8-6
编辑：Torsten 2023-8-6
Ok. "fsolve" might be able to cope with the problem, but usually there is something wrong with these tasks :-)
Torsten
Torsten 2023-8-6
编辑：Torsten 2023-8-6
@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
Yuting 2023-8-6
@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
Torsten 2023-8-6
编辑：Torsten 2023-8-6
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
Yuting 2023-8-6
编辑：Yuting 2023-8-6
@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
Torsten 2023-8-6
编辑：Torsten 2023-8-6
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
Walter Roberson 2023-8-7
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
Yuting 2023-8-7
@Torsten It would be something like this:
or using integral2() with xi as another variable needs to be integrated:
Yuting
Yuting 2023-8-7
@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))};

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Torsten
Torsten 2023-8-7

1 个投票

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

7 个评论
显示 5更早的评论 隐藏 5更早的评论

Yuting
Yuting 2023-8-7
Thanks!! The code is running without any error report now! Now my mission is to find the right initial guess values
Thanks again!
Torsten
Torsten 2023-8-7
... 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
Yuting 2023-8-8
Okay! I will definitely do it
Thanks again!
Best wishes,
Yuting
Yuting
Yuting 2023-8-10
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
Torsten 2023-8-10
编辑：Torsten 2023-8-10
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.
Torsten
Torsten 2023-8-11
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];

请先登录，再进行评论。

更多回答（1 个）

Walter Roberson
Walter Roberson 2023-8-6

1 个投票

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.

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Yuting
Yuting 2023-8-6
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
Walter Roberson 2023-8-6
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
Walter Roberson 2023-8-6
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
Yuting 2023-8-6
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...

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Conversion Between Symbolic and Numeric 的更多信息

产品

版本

R2022b

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by