Error: Limits of integration must be double or single scalars. While solving equations numerically with symbolic integral limits
33 views (last 30 days)
Show older comments
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:
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);
++++++++++++++++++++++++++
15 Comments
Accepted Answer
Torsten
on 7 Aug 2023
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 Comments
Torsten
on 11 Aug 2023
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];
More Answers (1)
Walter Roberson
on 6 Aug 2023
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 Comments
Walter Roberson
on 6 Aug 2023
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);
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!