Error : Not Enough Input Arguments at u1 = z(1)
Show older comments
Hello, I am new to matlab. I wrote this script by going through basics and by looking up answers in this community.
So, I am trying to solve six simultaneous non linear equations with six unknowns. I am using the fsolve method. I am working on parametric surface equations. I am getting the following error, Not Enough Input Arguments at u1 = z(1);
My Questions are:
1) How to solve this issue ?
2) How to get the values of the angles u1, v1, u2 and v2 in either radians or degrees?
function F=myFunction(z)
syms u1 v1 R r u2 v2 L t lambda Cx Cy Cz u0 v0 k df x ui vi;
%Six Unknowns
u1 = z(1);
v1 = z(2);
u2 = z(3);
v2 = z(4);
lambda = z(5);
t = z(6);
%Define Known Quantities
r=3; %Radius of the probe tip
R=5; %Radius of the Hemisphere
L=50; %Length of the Stylus Stem
k=1; %Number of the subprocess
df=0.001; %Deflection along the feed direction
%Define Surface 1 (Hemispherical Surface) in Parametric Form
S1=[R*cos(u1)*sin(v1);
R*sin(u1)*sin(v1);
R*cos(v1)];
%Define Surface 2(Probe tip) in Parametric Form
S2=[Cx+r*cos(u2)*sin(v2);
Cy+r*sin(u2)*sin(v2);
Cz+r*cos(v2)];
%Surface 1 Normal
n1 = cross(diff(S1,u1),diff(S1,v1));
%Normal Vector at Initial Point of Contact
n0 = subs(n1,[u1,v1],[deg2rad(0),deg2rad(50)]);
%Vector Vk in XY Plane
v0 = n0-[0;0;n0(3)];
%Initial Position of the Spindle
S0 = [(r+R)*cos(ui)*sin(vi);
(r+R)*sin(ui)*sin(vi);
(r+R)*cos(vi)+L];
%ui = deg2rad(ui);
%vi = deg2rad(vi);
S0 = subs(S0,[ui,vi],[deg2rad(0),deg2rad(50)]);
%Position of Spindle at any subprocess k
Sk = S0 + k*df*n0;
%Cx = Sk(1)+v0(1)*t;
%Cy = Sk(2)+v0(2)*t;
%Cz = Sk(3)-L;
S2=subs(S2,[Cx,Cy,Cz],[Sk(1)+v0(1)*t,Sk(2)+v0(2)*t,Sk(3)-L]);
%Surface 2 Normal
n2 = cross(diff(S2,u2),diff(S2,v2));
%Six Simultaneous Equations
F(1) = S1(1)-S2(1);
F(2) = S1(2)-S2(2);
F(3) = S1(3)-S2(3);
F(4) = n1(1)-lambda*n2(1);
F(5) = n1(2)-lambda*n2(2);
F(6) = n1(3)-lambda*n2(3);
F = [F(1),F(2),F(3),F(4),F(5),F(6)];
end
To run this function, I use the following code,
%guess values of six unknowns
zg = [0;0;0;0;1;1];
options = optimoptions('fmincon',...
'Algorithm','sqp','Display','iter','ConstraintTolerance',1e-12,'StepTolerance',1e-10);
z = fsolve(@myFunction, zg,options);
Thank you, If you need additional information or you think there is something missing please feel free to comment.
Accepted Answer
%guess values of six unknowns
zg = [0;0;0;0;1;1];
options = optimoptions('fsolve',...
'Algorithm','trust-region-dogleg','Display','iter','StepTolerance',1e-10);
z = fsolve(@myFunction, zg,options);
format long g
disp(z)
function F=myFunction(z)
syms u1 v1 R r u2 v2 L t lambda Cx Cy Cz u0 v0 k df x ui vi;
syms U1 U2 V1 V2
%Six Unknowns
u1 = z(1);
v1 = z(2);
u2 = z(3);
v2 = z(4);
lambda = z(5);
t = z(6);
%Define Known Quantities
r=3; %Radius of the probe tip
R=5; %Radius of the Hemisphere
L=50; %Length of the Stylus Stem
k=1; %Number of the subprocess
df=0.001; %Deflection along the feed direction
%Define Surface 1 (Hemispherical Surface) in Parametric Form
S1=[R*cos(U1)*sin(V1);
R*sin(U1)*sin(V1);
R*cos(V1)];
%Define Surface 2(Probe tip) in Parametric Form
S2=[Cx+r*cos(U2)*sin(V2);
Cy+r*sin(U2)*sin(V2);
Cz+r*cos(V2)];
%Surface 1 Normal
n1 = cross(diff(S1,U1),diff(S1,V1));
%Normal Vector at Initial Point of Contact
n0 = subs(n1,[U1,V1],[deg2rad(0),deg2rad(50)]);
%Vector Vk in XY Plane
v0 = n0-[0;0;n0(3)];
%Initial Position of the Spindle
S0 = [(r+R)*cos(ui)*sin(vi);
(r+R)*sin(ui)*sin(vi);
(r+R)*cos(vi)+L];
%ui = deg2rad(ui);
%vi = deg2rad(vi);
S0 = subs(S0,[ui,vi],[deg2rad(0),deg2rad(50)]);
%Position of Spindle at any subprocess k
Sk = S0 + k*df*n0;
%Cx = Sk(1)+v0(1)*t;
%Cy = Sk(2)+v0(2)*t;
%Cz = Sk(3)-L;
S2=subs(S2,[Cx,Cy,Cz],[Sk(1)+v0(1)*t,Sk(2)+v0(2)*t,Sk(3)-L]);
%Surface 2 Normal
N2 = cross(diff(S2,U2),diff(S2,V2));
n2 = subs(N2, [U2,V2], [u2,v2]);
%Six Simultaneous Equations
F(1) = S1(1)-S2(1);
F(2) = S1(2)-S2(2);
F(3) = S1(3)-S2(3);
F(4) = n1(1)-lambda*n2(1);
F(5) = n1(2)-lambda*n2(2);
F(6) = n1(3)-lambda*n2(3);
F = double(subs([F(1),F(2),F(3),F(4),F(5),F(6)], [U1,U2,V1,V2], [u1, u2, v1, v2]));
end
11 Comments
Muerus Rodrigues
on 20 Mar 2021
Walter Roberson
on 20 Mar 2021
In most cases you should not use symbolic expressions inside an fsolve objective function. You should instead call such a function once with symbolic arguements, and then use matlabFunction() to generate a numeric function for use with fsolve.
Your function involves, for example, taking the derivative of an expression with respect to a point. It is not useful to have to calculate the function and its derivative each time: you can be more efficient by building the function and taking its derivative once outside the objective function, and convert it to numeric anonymous function for use inside the objective.
There are some functions that are not defined except in Symbolic Toolbox, and there are some cases where you need to evaluate an expression to a lot of digits in order to get correct ratios, in which case you might need to duck into symbolic work for a few lines. It is, for example, common to take ratios of Bessel functions, but those explode to 10^1000 or 10^-1000 easily, so it just doesn't work to calculate them numerically, in which case Symbolic might be needed.
So it isn't "never" on using symbolic work within an fsolve objective function, but it is a case of "precompute what you can and turn that into numeric work."
Muerus Rodrigues
on 20 Mar 2021
Muerus Rodrigues
on 20 Mar 2021
Walter Roberson
on 20 Mar 2021
Edited: Walter Roberson
on 20 Mar 2021
When you use fsolve() it is not possible to put in constraints. However, you can use this adaptation.
You can change lb and ub to put in vectors of bounds.
syms Z [1 6]
myFunctionSym = myFunction(Z);
residueSym = sum(myFunctionSym.^2);
residue = matlabFunction(residueSym, 'vars', {Z.'});
%guess values of six unknowns
zg = [0;0;0;0;1;1];
options = optimoptions('fmincon',...
'Algorithm','sqp','Display','iter','ConstraintTolerance',1e-12,'StepTolerance',1e-10);
A = []; b = []; Aeq = []; beq = [];
lb = []; %vector of values of length 6 can be given here
ub = []; %vector of values of length 6 can be given here
nonlcon = [];
z = fmincon(residue, zg, A, b, Aeq, beq, lb, ub, nonlcon, options);
format long g
disp(z)
function F=myFunction(z)
syms u1 v1 R r u2 v2 L t lambda Cx Cy Cz u0 v0 k df x ui vi;
syms U1 U2 V1 V2
%Six Unknowns
u1 = z(1);
v1 = z(2);
u2 = z(3);
v2 = z(4);
lambda = z(5);
t = z(6);
%Define Known Quantities
r=3; %Radius of the probe tip
R=5; %Radius of the Hemisphere
L=50; %Length of the Stylus Stem
k=1; %Number of the subprocess
df=0.001; %Deflection along the feed direction
%Define Surface 1 (Hemispherical Surface) in Parametric Form
S1=[R*cos(U1)*sin(V1);
R*sin(U1)*sin(V1);
R*cos(V1)];
%Define Surface 2(Probe tip) in Parametric Form
S2=[Cx+r*cos(U2)*sin(V2);
Cy+r*sin(U2)*sin(V2);
Cz+r*cos(V2)];
%Surface 1 Normal
n1 = cross(diff(S1,U1),diff(S1,V1));
%Normal Vector at Initial Point of Contact
n0 = subs(n1,[U1,V1],[deg2rad(0),deg2rad(50)]);
%Vector Vk in XY Plane
v0 = n0-[0;0;n0(3)];
%Initial Position of the Spindle
S0 = [(r+R)*cos(ui)*sin(vi);
(r+R)*sin(ui)*sin(vi);
(r+R)*cos(vi)+L];
%ui = deg2rad(ui);
%vi = deg2rad(vi);
S0 = subs(S0,[ui,vi],[deg2rad(0),deg2rad(50)]);
%Position of Spindle at any subprocess k
Sk = S0 + k*df*n0;
%Cx = Sk(1)+v0(1)*t;
%Cy = Sk(2)+v0(2)*t;
%Cz = Sk(3)-L;
S2=subs(S2,[Cx,Cy,Cz],[Sk(1)+v0(1)*t,Sk(2)+v0(2)*t,Sk(3)-L]);
%Surface 2 Normal
N2 = cross(diff(S2,U2),diff(S2,V2));
n2 = subs(N2, [U2,V2], [u2,v2]);
%Six Simultaneous Equations
F(1) = S1(1)-S2(1);
F(2) = S1(2)-S2(2);
F(3) = S1(3)-S2(3);
F(4) = n1(1)-lambda*n2(1);
F(5) = n1(2)-lambda*n2(2);
F(6) = n1(3)-lambda*n2(3);
F = (subs([F(1),F(2),F(3),F(4),F(5),F(6)], [U1,U2,V1,V2], [u1, u2, v1, v2]));
end
Muerus Rodrigues
on 20 Mar 2021
Muerus Rodrigues
on 20 Mar 2021
Walter Roberson
on 20 Mar 2021
Notice that here I used the technique of calling the function once with symbolic arguments to get a symbolic expression, and then turning that into an anonymous handle with matlabFunction() .
Muerus Rodrigues
on 20 Mar 2021
Walter Roberson
on 20 Mar 2021
Ah, you would proceed a bit differently in that case. You would still create
myFunctionSym = myFunction(Z);
but you would then
myFunction_numeric = matlabFunction(myFunctionSym, 'vars', {Z.'});
nonlcon = @(z) deal(myFunction_numeric(z), []);
z = fmincon(Function_to_minimize, zg, A, b, Aeq, beq, lb, ub, nonlcon, options);
Muerus Rodrigues
on 20 Mar 2021
More Answers (1)
Muerus Rodrigues
on 24 Mar 2021
7 Comments
Walter Roberson
on 24 Mar 2021
How can you tell that the results are incorrect?
When I execute, I get
4.92017133636611e-12 3.51764035409061e-25 1.05521301865511e-32 2.26415762710886 -1.98782626363788e-24 1.00816899335998
which look reasonable. They are within the ub and lb, and the nonlinear equality constraints are at most 1e-17 . The function value is 4.55219293546628e-17 which seems pretty low.
rad2deg(z(1:4)) says that z(4) is 129.726676185693 degree, which is well within the limit of 180 degrees for it.
Muerus Rodrigues
on 25 Mar 2021
Walter Roberson
on 26 Mar 2021
There are multiple solutions. I used a different software package to find several exact solutions, but along the way I encounted bugs in that other package, so I have to go back and validate the solutions more carefully.
Two solutions that look legitimate are
Z1 = 0
Z2 = atan(1/1917*413437^(1/2)*3^(1/2))
Z3 = pi,
Z4 = -arctan(1/1917*413437^(1/2)*3^(1/2)) + pi
Z5 = -25/9
Z6 = -1/1000*413437^(1/2)+639/1000
Z1 = pi
Z2 = atan(1/1917*413437^(1/2)*3^(1/2))
Z3 = 0
Z4 = -arctan(1/1917*413437^(1/2)*3^(1/2)) + pi
Z5 = -25/9
Z6 = 1/1000*413437^(1/2)+639/1000
There is a third solution that looks self-consistent but which I need to cross-check.
The software package found another 7 solutions (10 total), some of which I might be able to modify into something workable, but others look just wrong so I am going to have to go back and work step-by-step to see what I can find.
Walter Roberson
on 26 Mar 2021
Edited: Walter Roberson
on 26 Mar 2021
Other solutions include
Z1 = anything within the bounds
Z2 = 0
Z3 = 0
Z4 = -atan((-1634563+1022400*3^(1/2))^(1/2)/(-800+639*3^(1/2))) + pi
Z5 = 0
Z6 = 1/1000*(-1634563+1022400*3^(1/2))^(1/2) + 639/1000
Z1 = anything within the bounds
Z2 = 0
Z3 = pi
Z4 = pi - acos(-5/3 + (213*sqrt(3))/160)
Z5 = 0
Z6 = -sqrt(-1634563 + 1022400*sqrt(3))/1000 + 639/1000
Walter Roberson
on 27 Mar 2021
All ten solution families that the other software package found, lead to an objective function value of exactly 0. But I have reason to believe that not all 10 solutions are consistent with the bounds.
Muerus Rodrigues
on 27 Mar 2021
Walter Roberson
on 27 Mar 2021
Edited: Walter Roberson
on 27 Mar 2021
Combining the above with some more results:
Pi = pi; arccos = @acos; arctan = @atan;
Z1 = 0, Z2 = 0, Z3 = 0, Z4 = Pi-arccos(-5/3+213/160*3^(1/2)), Z5 = 0, Z6 = 1/1000*(-1634563+1022400*3^(1/2))^(1/2)+639/1000
Z1 = 0, Z2 = 0, Z3 = Pi, Z4 = Pi-arccos(-5/3+213/160*3^(1/2)), Z5 = 0, Z6 = -1/1000*(-1634563+1022400*3^(1/2))^(1/2)+639/1000
Z1 = 0, Z2 = arctan(1/1917*413437^(1/2)*3^(1/2)), Z3 = pi, Z4 = -arctan(1/1917*413437^(1/2)*3^(1/2)) + pi, Z5 = -25/9, Z6 = -1/1000*413437^(1/2)+639/1000
Z1 = Pi, Z2 = 0, Z3 = 0, Z4 = Pi-arccos(-5/3+213/160*3^(1/2)), Z5 = 0, Z6 = 1/1000*(-1634563+1022400*3^(1/2))^(1/2)+639/1000
Z1 = Pi, Z2 = 0, Z3 = Pi, Z4 = Pi-arccos(-5/3+213/160*3^(1/2)), Z5 = 0, Z6 = -1/1000*(-1634563+1022400*3^(1/2))^(1/2)+639/1000
Z1 = pi, Z2 = arctan(1/1917*413437^(1/2)*3^(1/2)), Z3 = 0, Z4 = -arctan(1/1917*413437^(1/2)*3^(1/2)) + pi, Z5 = -25/9, Z6 = 1/1000*413437^(1/2)+639/1000
Z1 = Z1, Z2 = 0, Z3 = 0, Z4 = Pi-arccos(-5/3+213/160*3^(1/2)), Z5 = 0, Z6 = 1/1000*(-1634563+1022400*3^(1/2))^(1/2)+639/1000
Z1 = Z1, Z2 = 0, Z3 = 0, Z4 = -arctan((-1634563+1022400*3^(1/2))^(1/2)/(-800+639*3^(1/2))) + pi, Z5 = 0, Z6 = 1/1000*(-1634563+1022400*3^(1/2))^(1/2) + 639/1000
Z1 = Z1, Z2 = 0, Z3 = Pi, Z4= Pi-arccos(-5/3+213/160*3^(1/2)), Z5 = 0, Z6 = -1/1000*(-1634563+1022400*3^(1/2))^(1/2)+639/1000
Z1 = Z1, Z2 = 0, Z3 = pi, Z4 = pi - arccos(-5/3 + (213*sqrt(3))/160), Z5 = 0, Z6 = -sqrt(-1634563 + 1022400*sqrt(3))/1000 + 639/1000
All of these lead to an f value of exactly 0, which is the minimum possible for that function.
The notation Z1 = Z1 indicates that Z1 can be any value within its permitted range.
Categories
Find more on Collect Coverage for Models in Help Center and File Exchange
on 20 Mar 2021
on 27 Mar 2021
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
