MATLAB: Problem with fsolve for non linear system

fsolveMATLABnonlinear system

Hi there, I'm trying to solve this system with fsolver:
F1 = P(1) - P(3) - (- dV1(i+1)*vel(i+1)/D - sqrt(abs(P(1) - P(2)))*sign(P(1) - P(2)))^2 * sign(- dV1(i+1)*vel(i+1)/D - sqrt(abs(P(1) - P(2)))*sign(P(1) - P(2)));
F2 = P(1) - P(2) - (dV2(i+1)*vel(i+1)/D + sqrt(abs(P(2) - P(3)))*sign(P(2) - P(3)))^2 * sign(dV2(i+1)*vel(i+1)/D + sqrt(abs(P(2) - P(3)))*sign(P(2) - P(3)));
F3 = P(1) + ((V1(i+1)*(P(1) - P1(i,1)) + V2(i+1)*(P(2) - P2(i,1)) + V3(i+1)*(P(3) - P3(i,1)))/Dtheta(i+1) + P(2)*dV2(i+1) + P(3)*dV3(i+1))/dV1(i+1);
F = [F1 ; F2 ; F3];
P0 = [P1(i,1); P2(i,1); P3(i,1)];
options = optimoptions('fsolve','Display','iter');
[P,fval] = fsolve(F,P0,options)
And it returns:
    Error using lsqfcnchk (line 108)
  FUN must be a function, a valid string expression, or an inline function object.
Error in fsolve (line 197)
    funfcn = lsqfcnchk(FUN,'fsolve',length(varargin),funValCheck,gradflag);
Error in wankinc (line 81)
    [P,fval] = fsolve(F,P0,options)
Could anyone tell me what the problem is? I don't see why the values of F1, F2 and F3 are not valid.
Thank you!
PS: all the variables within the functions are calculated in the same code, previously. The complete code is copied below, if anyone finds it necessary:
e = 3/1000;     %geometry
R = 21/1000;
W = 5/1000;
G=3*sqrt(3)*2000*W*R^2*(3*R^2/8 + e^2)/4;   %Rotational inertia
P1(1,1) = 1.01325e5 - 1e3;    %starting fluid pressures, from equilibrium
P2(1,1) = 1.01325e5 + 1e3;
P3(1,1) = 1.01325e5;
P= [P1(1,1) P2(1,1) P3(1,1)];
dP1(1) = 0; %pressure change over one timestep
dP2(1) = 0;
dP3(1) = 0;
X(1) = 0;   %pressure gradient between chambers. Properly defined ahead.
Y(1) = 0;
Z(1) = 0;
T_app(1) = 1;   %T_app(1) is the torque to be ultimately imposed on the system
T_app(2) = T_app(1)/20;   %T_app(2) is a variable torque that increases from 0 up to the imposed torque
T_res(1) = 0;   %Torque caused by the pressure in the chambers
T = zeros(1010);
theta(1) = 0;   %Displacement over the entire simultaion
Dtheta(1) = 0;  %Displacement over one time step
vel(1) = 0;     %Angular velocity
alpha(1) = 0;   %Angular acceleration
A = pi*(e^2 + R^2/3) - sqrt(3)*R^2/4;       %A, B, C and D are constants defined here for computational speed
B = 3*sqrt(3)*e*R/4;
C = 3*sqrt(3)*e*R/2;
D = 0.6 * pi*0.1^2/4 * sqrt(2/1000);      %This is the flow constant for the orifices
E = 1e3;                                      %Error of 1kPa is accepted before continuing to the next step
f = zeros(1500,4);
f(:,:) = 2*E;
V1(1) = A - B * (sqrt(3)*sin(2*(theta(1) - 2*pi/3) + cos(2*(theta(1) - 2*pi/3))));      %Chamber volume, a function of theta
dV1(1) = C*(sin(2*(theta(1) - 2*pi/3)) - sqrt(3)*cos(2*(theta(1) - 2*pi/3)));           %dV/dtheta
V2(1) = A - B * (sqrt(3)*sin(2*theta(1)) + cos(2*theta(1)));
dV2(1) = C*(sin(2*theta(1)) - sqrt(3)*cos(2*theta(1)));
V3(1) = A - B * (sqrt(3)*sin(2*(theta(1) + 2*pi/3) + cos(2*(theta(1) + 2*pi/3))));
dV3(1) = C*(sin(2*(theta(1) + 2*pi/3)) - sqrt(3)*cos(2*(theta(1) + 2*pi/3)));
ts = 1e-3;                                  %Time step of 0.1 ms
alpha(1) = T_app(2)/G;
for i = 1:1:1
      if T_app(2)<T_app(1)                         %Ramps up the torque in each iteration
          T_app(2) = T_app(2) + T_app(1)/20;
      end
      vel(i+1) = vel(i) + alpha(i)*ts;
      theta(i+1) = theta(i) + vel(i)*ts + alpha(i)*ts^2/2;
      Dtheta(i+1) = theta(i+1) - theta(i);
      V1(i+1) = A - B * (sqrt(3)*sin(2*(theta(i+1)) - 2*pi/3) + cos(2*(theta(i+1) - 2*pi/3)));
      V2(i+1) = A - B * (sqrt(3)*sin(2*theta(i+1)) + cos(2*theta(i+1)));
      V3(i+1) = A - B * (sqrt(3)*sin(2*(theta(i+1)) + 2*pi/3) + cos(2*(theta(i+1) + 2*pi/3)));
      dV1(i+1) = C*(sin(2*(theta(i+1) - 2*pi/3)) - sqrt(3)*cos(2*(theta(i+1) - 2*pi/3)));
      dV2(i+1) = C*(sin(2*theta(i+1)) - sqrt(3)*cos(2*theta(i+1)));
      dV3(i+1) = C*(sin(2*(theta(i+1) + 2*pi/3)) - sqrt(3)*cos(2*(theta(i+1) + 2*pi/3)));
      F1 = P(1) - P(3) - (- dV1(i+1)*vel(i+1)/D - sqrt(abs(P(1) - P(2)))*sign(P(1) - P(2)))^2 * sign(- dV1(i+1)*vel(i+1)/D - sqrt(abs(P(1) - P(2)))*sign(P(1) - P(2)));
      F2 = P(1) - P(2) - (dV2(i+1)*vel(i+1)/D + sqrt(abs(P(2) - P(3)))*sign(P(2) - P(3)))^2 * sign(dV2(i+1)*vel(i+1)/D + sqrt(abs(P(2) - P(3)))*sign(P(2) - P(3)));
      F3 = P(1) + ((V1(i+1)*(P(1) - P1(i,1)) + V2(i+1)*(P(2) - P2(i,1)) + V3(i+1)*(P(3) - P3(i,1)))/Dtheta(i+1) + P(2)*dV2(i+1) + P(3)*dV3(i+1))/dV1(i+1);
      F = [F1 ; F2 ; F3];
      P0 = [P1(i,1); P2(i,1); P3(i,1)];
      options = optimoptions('fsolve','Display','iter');
      [P,fval] = fsolve(F,P0,options)
      P1(i+1,1) = P(1);
      P2(i+1,1) = P(2);
      P3(i+1,1) = P(3);
      T_res(i+1) = sqrt(3)*R*W*e * (-P1(i+1,1)*cos(2*theta(i+1) - pi/6) + P2(i+1,1)*cos(2*theta(i+1) + pi/6) + P3(i+1,1)*sin(2*theta(i+1)));
      T(i+1) = T_app(2) - T_res(i+1);
      alpha(i+1) = T(i+1)/G;
end

Best Answer

Do this:
- in your main function, after calculating P1,P2,P3, put:
P0 = [P1(i,1); P2(i,1); P3(i,1)];
options = optimoptions('fsolve','Display','iter');
[P,fval] = fsolve(@myfun,P0,options)
- create a second function
function F = myfun(P)
F1 = P(1) - P(3) - (- dV1(i+1)*vel(i+1)/D - sqrt(abs(P(1) - P(2)))*sign(P(1) - P(2)))^2 * sign(- dV1(i+1)*vel(i+1)/D - sqrt(abs(P(1) - P(2)))*sign(P(1) - P(2)));
F2 = P(1) - P(2) - (dV2(i+1)*vel(i+1)/D + sqrt(abs(P(2) - P(3)))*sign(P(2) - P(3)))^2 * sign(dV2(i+1)*vel(i+1)/D + sqrt(abs(P(2) - P(3)))*sign(P(2) - P(3)));
F3 = P(1) + ((V1(i+1)*(P(1) - P1(i,1)) + V2(i+1)*(P(2) - P2(i,1)) + V3(i+1)*(P(3) - P3(i,1)))/Dtheta(i+1) + P(2)*dV2(i+1) + P(3)*dV3(i+1))/dV1(i+1);
F = [F1 ; F2 ; F3];
To pass the parameters, you can either use a global or do:
[P,fval] = fsolve(@(x)myfun(x,param1,param2),P0,options)
function F = myfun(P,param1,param2)
Related Question