Problem with fsolve for non linear 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

Answer 1

[Sara]

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)

[Julio Caye]

when you tell me to create a second function, is it really creating a new function file and saving it as a ".m" file, or to insert those lines in the main body?

[Sara]

It depends on how you write your code. If the main body is in a function itself, i.e. the .m file starts with the word function filename, then you can add the new function at the bottom. Otherwise, you'll have to create a new m file and put it in the same folder as the file calling such function

[Julio Caye]

thanks a lot! The results are a bit non physical, but that's most likely my mistake while deriving the equations... hahaha

Thanks again!

[Sara]

There I can't help :)

[Julio Caye]

Can I ask just one more thing, please? How do I make it not plot out the number of iterations and error and stuff? I have to run that system some 10 thousand steps and that information kinda overwhelms me haha

[Julio Caye]

nevermind, found the option to do it here (and btw, the eqations were right, just had to reduce the timestep a bit)

thanks again!