MATLAB: FSOLVE requires all values returned by functions to be of data type double.

fsolve doubleOptimization ToolboxSymbolic Math Toolboxsymbolic vs numeric

dear all i write code like below that produce 9*1 matrix that have 9 sym variable x1,x2,………,x9 i m using fsolve to calculate these variable :
      function F = consnt4(x)
      m=1.669e-3;
      jx=5.743e-14;
      E=2.1e11;
      O=12.5e-3;
      jz=1.149e-13;
      G=8e10;
      L=3e-3;
      g=1;%9.81;
      N=3;
      T1=5;
      T2=0;
      T3=0;
      syms  real
      q    = sym('q',    [1 N+1], 'real');
      eb   = sym('eb',   [1 N+1], 'real');
      teta = sym('teta', [1 N+1], 'real');
      tetad = sym('tetad', [1 N+1], 'real');
      ebd   = sym('ebd',   [1 N+1], 'real');
      qd    = sym('qd',    [1 N+1], 'real');
      x      =sym('x',    [1 3*N+1], 'real');
    r(:,:,1)=[1*O;0;0;];
    r(:,:,2)=[(-1/2)*O;(sqrt(3))*O/2;0];
    r(:,:,3)=[(-1/2)*O;-(sqrt(3))*O/2;0];
    for i=2:N+1
        pl(:,:,i)=L*[(cos(q(i))*(1-cos(teta(i))))/teta(i);(sin(q(i))*(1-cos(teta(i))))/teta(i) ; sin(teta(i))/teta((i))];
        Rl(:,:,i)=[cos(q(i))*cos(teta(i))*cos(eb(i)-q(i))-sin(q(i))*sin(eb(i)-q(i)) -cos(q(i))*cos(teta(i))*sin(eb(i)-q(i))-sin(q(i))*cos(eb(i)-q(i)) cos(q(i))*sin(teta(i)); sin(q(i))*cos(teta(i))*cos(eb(i)-q(i))+cos(q(i))*sin(eb(i)-q(i)) -sin(q(i))*cos(teta(i))*sin(eb(i)-q(i))+cos(q(i))*cos(eb(i)-q(i)) sin(q(i))*sin(teta(i));-sin(teta(i))*cos(eb(i)-q(i)) sin(teta(i))*sin(eb(i)-q(i)) cos(teta(i))];
        R(:,:,2)=Rl(:,:,2);
        R(:,:,1)=[1 0 0;0 1 0;0 0 1];
        p(:,:,2)=pl(:,:,2);
    end  
    for i=3:N+1
        R(:,:,i)=R(:,:,i-1)*Rl(:,:,i);
        p(:,:,i)=p(:,:,i-1)+R(:,:,i-1)*pl(:,:,i);
    end
    for i=2:N+1
        wl(:,:,i)=[-sin(q(i))*tetad(i)-cos(q(i))*sin(teta(i))*cos(teta(i))*qd(i)+cos(q(i))*sin(teta(i))*ebd(i);cos(q(i))*tetad(i)-sin(q(i))*sin(teta(i))*cos(teta(i))*qd(i)+sin(q(i))*sin(teta(i))*ebd(i);((sin(teta(i)))^2)*qd(i)+cos(teta(i))*ebd(i)];
        w(:,:,2)=wl(:,:,2);
    end
    for i=3:N+1
        w(:,:,i)= w(:,:,i-1)+R(:,:,i-1)*wl(:,:,i);
    end
    for i=2:N+1
        for j=2:i
           wc(:,:,i)=collect(w(:,:,i),[qd(2),tetad(2),ebd(2),qd(3),tetad(3),ebd(3),qd(4),tetad(4),ebd(4)]) ;  %   ;,qd(5),tetad(5),ebd(5),qd(6),tetad(6),ebd(6),qd(7),tetad(7),ebd(7),qd(8),tetad(8),ebd(8),qd(9),tetad(9),ebd(9)]);
        end
    end
    for i=2:N+1
        for j=1:size(w(:,:,i),1)
            [wtemp,vtemp]=coeffs(wc(j,:,i),[qd(2),tetad(2),ebd(2),qd(3),tetad(3),ebd(3),qd(4),tetad(4),ebd(4)]);%  ,qd(5),tetad(5),ebd(5),qd(6),tetad(6),ebd(6),qd(7),tetad(7),ebd(7),qd(8),tetad(8),ebd(8),qd(9),tetad(9),ebd(9)]);

            [~,idx]=ismember(vtemp,[qd(2),tetad(2),ebd(2),qd(3),tetad(3),ebd(3),qd(4),tetad(4),ebd(4)]);%  ,qd(5),tetad(5),ebd(5),qd(6),tetad(6),ebd(6),qd(7),tetad(7),ebd(7),qd(8),tetad(8),ebd(8),qd(9),tetad(9),ebd(9)]);
            wk(j,idx,i)=wtemp;
        end
    end
    for i=2:N+1
      pd(:,:,i)=[teta(i)*((-qd((i))*sin(q(i))*(1-cos(teta(i)))+tetad(i)*sin(teta(i))*cos(q(i)))+tetad(i)*(L*cos(q(i))*(1-cos(teta(i)))))/(teta((i))^2); teta(i)*((qd((i))*cos(q(i))*(1-cos(teta(i)))+tetad(i)*sin(teta(i))*sin(q(i)))+tetad(i)*(L*sin(q(i))*(1-cos(teta(i)))))/(teta((i))^2);  (teta(i)*tetad(i)*cos(teta(i))-tetad(i)*L*sin(teta(i)))/(teta(i)^2)];
      v(:,:,2)=pd(:,:,2);
     end
     for i=3:N+1
         v(:,:,i)= v(:,:,i-1)+cross(w(:,:,i-1),R(:,:,i-1)*pl(:,:,i))+R(:,:,i-1)*pd(:,:,i);
     end
    for i=2:N+1
        for j=2:i
            vc(:,:,i)=collect(v(:,:,i),[qd(2),tetad(2),ebd(2),qd(3),tetad(3),ebd(3),qd(4),tetad(4),ebd(4)]);%,qd(5),tetad(5),ebd(5),qd(6),tetad(6),ebd(6),qd(7),tetad(7),ebd(7),qd(8),tetad(8),ebd(8),qd(9),tetad(9),ebd(9)]);


        end
    end
    for i=2:N+1
        for j=1:size(v(:,:,i),1)
            [vwtemp,vvtemp]=coeffs(vc(j,:,i),[qd(2),tetad(2),ebd(2),qd(3),tetad(3),ebd(3),qd(4),tetad(4),ebd(4)]);%,qd(5),tetad(5),ebd(5),qd(6),tetad(6),ebd(6),qd(7),tetad(7),ebd(7),qd(8),tetad(8),ebd(8),qd(9),tetad(9),ebd(9)]);
            [~,idx]=ismember(vvtemp,[qd(2),tetad(2),ebd(2),qd(3),tetad(3),ebd(3),qd(4),tetad(4),ebd(4)]);%,qd(5),tetad(5),ebd(5),qd(6),tetad(6),ebd(6),qd(7),tetad(7),ebd(7),qd(8),tetad(8),ebd(8),qd(9),tetad(9),ebd(9)]);
            vk(j,idx,i)=vwtemp;  
        end
                vk(:,3*(i-1),i)=[0;0;0];
    end 
    for i=N+1
      Mb(:,:,i)=E*jx*((teta(i))/L)*R(:,:,i-1)*[-sin(q(i));cos(q(i));0];
      Mt(:,:,i)=-G*jz*eb(i)*R(:,1,i)/L;
      Me(:,:,i)=Mt(:,:,i)-Mb(:,:,i);
    end
    for i=2:N
    Mb(:,:,i)=E*jx*((teta(i))/L)*R(:,:,i-1)*[-sin(q(i));cos(q(i));0];
    end
    for i=2:N
    Mt(:,:,i)=-G*jz*eb(i)*R(:,1,i)/L;
    end
    for i=3:N+1
    Mtt(:,:,i)=G*jz*eb(i)*R(:,1,i)/L;
    end
    for i=2:N
        MT(:,:,i)=Mt(:,:,i)+Mtt(:,:,i+1);
    end
    for  i=2:N
        Me(:,:,i)=Mb(:,:,i+1)-Mb(:,:,i)+MT(:,:,i);
    end
    for i=2:N+1
    Fg(:,:,i)=-m*g*[1;0;0];
    end
    for i=2
    ph(:,1,i)=[pl(:,:,i)+R(:,:,i)*r(:,:,1)-r(:,:,1)];
    ph(:,2,i)=[pl(:,:,i)+R(:,:,i)*r(:,:,2)-r(:,:,2)];
    ph(:,3,i)=[pl(:,:,i)+R(:,:,i)*r(:,:,3)-r(:,:,3)];
    end
    for i=3:N+1
    ph(:,1,i)=[R(:,:,i-1)*pl(:,:,i)+R(:,:,i)*r(:,:,1)-R(:,:,i-1)*r(:,:,1)];
    ph(:,2,i)=[R(:,:,i-1)*pl(:,:,i)+R(:,:,i)*r(:,:,2)-R(:,:,i-1)*r(:,:,2)];
    ph(:,3,i)=[R(:,:,i-1)*pl(:,:,i)+R(:,:,i)*r(:,:,3)-R(:,:,i-1)*r(:,:,3)];
    end
    for i=2:N+1
    ap(:,:,i)=[sqrt(ph(1,1,i)^2+ph(2,1,i)^2+ph(3,1,i)^2);sqrt(ph(1,2,i)^2+ph(2,2,i)^2+ph(3,2,i)^2);sqrt(ph(1,3,i)^2+ph(2,3,i)^2+ph(3,3,i)^2)];
    f(:,:,i)=[ph(:,1,i)/ap(1,1,i) ph(:,2,i)/ap(2,1,i) ph(:,3,i)/ap(3,1,i)];
    end
    for i=N+1
    Fc(:,:,i)=[-T1*f(:,1,i) -T2*f(:,2,i) -T3*f(:,3,i)];
    end
    for i=2:N
        Fc(:,:,i)=[T1*f(:,1,i+1) T2*f(:,2,i+1) T3*f(:,3,i+1)]+[-T1*f(:,1,i) -T2*f(:,2,i) -T3*f(:,3,i)];      
    end
    for i=2:N+1
       fa(:,:,i)=Fc(:,1,i)+Fc(:,2,i)+Fc(:,3,i);
    end
    for i=2:N+1
        Ma(:,:,i)=cross(r(:,:,1),Fc(:,1,i))+cross(r(:,:,2),Fc(:,2,i))+cross(r(:,:,3),Fc(:,3,i));
    end
    for i=2:N+1
    feq(:,:,i)=fa(:,:,i); +Fg(:,:,i);
    Meq(:,:,i)=Ma(:,:,i); + Me(:,:,i);
    end
    for i=2:N+1
        for j=2:N+1
       vk(:,:,i)=subs(vk(:,:,i),q(j),x(j-1));
       wk(:,:,i)=subs(wk(:,:,i),q(j),x(j-1));
       Meq(:,:,i)=subs(Meq(:,:,i),q(j),x(j-1));
       feq(:,:,i)=subs(feq(:,:,i),q(j),x(j-1));
        end
    end
    for i=2:N+1
        for j=2:N+1
       vk(:,:,i)=subs(vk(:,:,i),eb(j),x(j+2));
       wk(:,:,i)=subs(wk(:,:,i),eb(j),x(j+2));
       Meq(:,:,i)=subs(Meq(:,:,i),eb(j),x(j+2));
       feq(:,:,i)=subs(feq(:,:,i),eb(j),x(j+2));
        end
    end
    for i=2:N+1
        for j=2:N+1
       vk(:,:,i)=subs(vk(:,:,i),teta(j),x(j+5));
       wk(:,:,i)=subs(wk(:,:,i),teta(j),x(j+5));
       Meq(:,:,i)=subs(Meq(:,:,i),teta(j),x(j+5));
       feq(:,:,i)=subs(feq(:,:,i),teta(j),x(j+5));
        end
    end
    for j=1:(3*N)
        for i=2:N+1
         F(j,:,:)= dot(Meq(:,:,i),wk(:,j,i))+dot(feq(:,:,i),vk(:,j,i));
        end
    end
and after that i try
x0=[0;0];
options = optimoptions('fsolve','Display','iter');
[x,fval] = fsolve(@myfunn,x0,options)
but i got error : FSOLVE requires all values returned by functions to be of data type double. cod that you see seems to be complex but it is just some simple multiplication of matrix and there is no integration i really appreciated if you could help it is the matter of life and death

Best Answer

If you want to solve a system of equations symbolically or using symbolic variables, use the solve function instead of fsolve. The fsolve function requires that the output of the function you specify as the first input be double precision numbers, not symbolic expressions.
If you must use fsolve with symbolic variables inside the function, you must convert those symbolic expressions into numbers using double before returning the numbers (not the symbolic expressions) from the function.
As a very simple example, put all this code into a function example373677.m and run it:
function y = example373677
y = fsolve(@mysubfun, 3);
 
function z = mysubfun(q)
x = sym('x');
r = x.^3-2*x-5;
z = double(subs(r, q));
Of course, with solve this is easier.
syms x
y = solve(x.^3-2*x-5)
You may need to use vpa to approximate y numerically.
vpa(y)
One of the solutions returned in y should match the result you received from running the function.
Related Question