MATLAB: As a new guy, I can not find the problem for the following code. It would be greatly appriciated if you guys can give out the answer. The specific error and code are showed in the body part.

??? The following error occurred converting from sym to double:
Error using ==> sym.double at 29
DOUBLE cannot convert the input expression into a double array.
If the input expression contains a symbolic variable, use the VPA function instead.
Error in ==> method2 at 244
    J(j,i)=diff(f(j),var(i));
function var1=method2()
syms x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24 x25 x26 x27 x28 x29 x30 x31 x32 x33 x34 x35 x36 x37 x38 x39 x40 real
var=[x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21 x22 x23 x24 x25 x26 x27 x28 x29 x30 x31 x32 x33 x34 x35 x36 x37 x38 x39 x40];
lamda_o=2;
lamda_i=0;
Fa=1000;
Fr=0;
M=7.83e3*4/3*3.1415926*(10.7e-3/2)^3;
vb=0.3;
vi=0.3;
vo=0.3;
Eb=219e9;
Ei=219e9;
Eo=219e9;
ri=0.515*10.7e-3;
ro=0.525*10.7e-3;
D=10.7e-3;
w=2*3.1415926*1000/60;
r=10.7/77.5;
Dm=77.5e-3;
J=0.4*M*(D/2)^2;
fo=0.525;
fi=0.515;
a0=0.44;
Eii=1/((1-vb^2)/Eb+(1-vi^2)/Ei);
Eoo=1/((1-vb^2)/Eb+(1-vo^2)/Eo);
BD=(fo+fi-1)*D;
fai_1=0;
fai_2=0.3491;
fai_3=0.6981;
fai_4=1.0472;
fai_5=1.3963;
fai_6=1.7453;
fai_7=2.0944;
fai_8=2.4435;
fai_9=2.7925;
fai_10=3.1416;
fai_11=3.4907;
fai_12=3.8397;
fai_13=4.1888;
fai_14=4.5379;
fai_15=4.8869;
fai_16=5.2360;
fai_17=5.5851;
fai_18=5.9341;
fai_19=6.2832;
sigemao1=((BD*cos(a0)+x39*cos(fai_1))*sin(x1)-((fo-0.5)*D)*cos(x20)*sin(x1)-(BD*sin(a0)+x40)*cos(x1)+((fo-0.5)*D)*sin(x20)*cos(x1))/...
(sin(x20)*cos(x1)-cos(x20)*sin(x1));
sigemao2=((BD*cos(a0)+x39*cos(fai_2))*sin(x2)-((fo-0.5)*D)*cos(x21)*sin(x2)-(BD*sin(a0)+x40)*cos(x2)+((fo-0.5)*D)*sin(x21)*cos(x2))/...
(sin(x21)*cos(x2)-cos(x21)*sin(x2));
sigemao3=((BD*cos(a0)+x39*cos(fai_3))*sin(x3)-((fo-0.5)*D)*cos(x22)*sin(x3)-(BD*sin(a0)+x40)*cos(x3)+((fo-0.5)*D)*sin(x22)*cos(x3))/...
(sin(x22)*cos(x3)-cos(x22)*sin(x3));
sigemao4=((BD*cos(a0)+x39*cos(fai_4))*sin(x4)-((fo-0.5)*D)*cos(x23)*sin(x4)-(BD*sin(a0)+x40)*cos(x4)+((fo-0.5)*D)*sin(x23)*cos(x4))/...
(sin(x23)*cos(x4)-cos(x23)*sin(x4));
sigemao5=((BD*cos(a0)+x39*cos(fai_5))*sin(x5)-((fo-0.5)*D)*cos(x24)*sin(x5)-(BD*sin(a0)+x40)*cos(x5)+((fo-0.5)*D)*sin(x24)*cos(x5))/...
(sin(x24)*cos(x5)-cos(x24)*sin(x5));
sigemao6=((BD*cos(a0)+x39*cos(fai_6))*sin(x6)-((fo-0.5)*D)*cos(x25)*sin(x6)-(BD*sin(a0)+x40)*cos(x6)+((fo-0.5)*D)*sin(x25)*cos(x6))/...
(sin(x25)*cos(x6)-cos(x25)*sin(x6));
sigemao7=((BD*cos(a0)+x39*cos(fai_7))*sin(x7)-((fo-0.5)*D)*cos(x26)*sin(x7)-(BD*sin(a0)+x40)*cos(x7)+((fo-0.5)*D)*sin(x26)*cos(x7))/...
(sin(x26)*cos(x7)-cos(x26)*sin(x7));
sigemao8=((BD*cos(a0)+x39*cos(fai_8))*sin(x8)-((fo-0.5)*D)*cos(x27)*sin(x8)-(BD*sin(a0)+x40)*cos(x8)+((fo-0.5)*D)*sin(x27)*cos(x8))/...
(sin(x27)*cos(x8)-cos(x27)*sin(x8));
sigemao9=((BD*cos(a0)+x39*cos(fai_9))*sin(x9)-((fo-0.5)*D)*cos(x28)*sin(x9)-(BD*sin(a0)+x40)*cos(x9)+((fo-0.5)*D)*sin(x28)*cos(x9))/...
(sin(x28)*cos(x9)-cos(x28)*sin(x9));
sigemao10=((BD*cos(a0)+x39*cos(fai_10))*sin(x10)-((fo-0.5)*D)*cos(x29)*sin(x10)-(BD*sin(a0)+x40)*cos(x10)+((fo-0.5)*D)*sin(x29)*cos(x10))/...
(sin(x29)*cos(x10)-cos(x29)*sin(x10));
sigemao11=((BD*cos(a0)+x39*cos(fai_11))*sin(x11)-((fo-0.5)*D)*cos(x30)*sin(x11)-(BD*sin(a0)+x40)*cos(x11)+((fo-0.5)*D)*sin(x30)*cos(x11))/...
(sin(x30)*cos(x11)-cos(x30)*sin(x11));
sigemao12=((BD*cos(a0)+x39*cos(fai_12))*sin(x12)-((fo-0.5)*D)*cos(x31)*sin(x12)-(BD*sin(a0)+x40)*cos(x12)+((fo-0.5)*D)*sin(x31)*cos(x12))/...
(sin(x31)*cos(x12)-cos(x31)*sin(x12));
sigemao13=((BD*cos(a0)+x39*cos(fai_13))*sin(x13)-((fo-0.5)*D)*cos(x32)*sin(x13)-(BD*sin(a0)+x40)*cos(x13)+((fo-0.5)*D)*sin(x32)*cos(x13))/...
(sin(x32)*cos(x13)-cos(x32)*sin(x13));
sigemao14=((BD*cos(a0)+x39*cos(fai_14))*sin(x14)-((fo-0.5)*D)*cos(x33)*sin(x14)-(BD*sin(a0)+x40)*cos(x14)+((fo-0.5)*D)*sin(x33)*cos(x14))/...
(sin(x33)*cos(x14)-cos(x33)*sin(x14));
sigemao15=((BD*cos(a0)+x39*cos(fai_15))*sin(x15)-((fo-0.5)*D)*cos(x34)*sin(x15)-(BD*sin(a0)+x40)*cos(x15)+((fo-0.5)*D)*sin(x34)*cos(x15))/...
(sin(x34)*cos(x15)-cos(x34)*sin(x15));
sigemao16=((BD*cos(a0)+x39*cos(fai_16))*sin(x16)-((fo-0.5)*D)*cos(x35)*sin(x16)-(BD*sin(a0)+x40)*cos(x16)+((fo-0.5)*D)*sin(x35)*cos(x16))/...
(sin(x35)*cos(x16)-cos(x35)*sin(x16));
sigemao17=((BD*cos(a0)+x39*cos(fai_17))*sin(x17)-((fo-0.5)*D)*cos(x36)*sin(x17)-(BD*sin(a0)+x40)*cos(x17)+((fo-0.5)*D)*sin(x36)*cos(x17))/...
(sin(x36)*cos(x17)-cos(x36)*sin(x17));
sigemao18=((BD*cos(a0)+x39*cos(fai_18))*sin(x18)-((fo-0.5)*D)*cos(x37)*sin(x18)-(BD*sin(a0)+x40)*cos(x18)+((fo-0.5)*D)*sin(x37)*cos(x18))/...
(sin(x37)*cos(x18)-cos(x37)*sin(x18));
sigemao19=((BD*cos(a0)+x39*cos(fai_19))*sin(x19)-((fo-0.5)*D)*cos(x38)*sin(x19)-(BD*sin(a0)+x40)*cos(x19)+((fo-0.5)*D)*sin(x38)*cos(x19))/...
(sin(x38)*cos(x19)-cos(x38)*sin(x19));
sigemai1=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao1)*sin(x20))/sin(x1)-(fi-0.5)*D;
sigemai2=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao2)*sin(x21))/sin(x2)-(fi-0.5)*D;
sigemai3=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao3)*sin(x22))/sin(x3)-(fi-0.5)*D;
sigemai4=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao4)*sin(x23))/sin(x4)-(fi-0.5)*D;
sigemai5=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao5)*sin(x24))/sin(x5)-(fi-0.5)*D;
sigemai6=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao6)*sin(x25))/sin(x6)-(fi-0.5)*D;
sigemai7=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao7)*sin(x26))/sin(x7)-(fi-0.5)*D;
sigemai8=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao8)*sin(x27))/sin(x8)-(fi-0.5)*D;
sigemai9=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao9)*sin(x28))/sin(x9)-(fi-0.5)*D;
sigemai10=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao10)*sin(x29))/sin(x10)-(fi-0.5)*D;
sigemai11=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao11)*sin(x30))/sin(x11)-(fi-0.5)*D;
sigemai12=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao12)*sin(x31))/sin(x12)-(fi-0.5)*D;
sigemai13=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao13)*sin(x32))/sin(x13)-(fi-0.5)*D;
sigemai14=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao14)*sin(x33))/sin(x14)-(fi-0.5)*D;
sigemai15=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao15)*sin(x34))/sin(x15)-(fi-0.5)*D;
sigemai16=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao16)*sin(x35))/sin(x16)-(fi-0.5)*D;
sigemai17=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao17)*sin(x36))/sin(x17)-(fi-0.5)*D;
sigemai18=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao18)*sin(x37))/sin(x18)-(fi-0.5)*D;
sigemai19=((BD*sin(a0))+x40-((fo-0.5)*D+sigemao19)*sin(x38))/sin(x19)-(fi-0.5)*D;
Mg1=J/r*((1-r*cos(x1))/(1+cos(x1-x20)))^2*w^2*sin(x20);
Mg2=J/r*((1-r*cos(x2))/(1+cos(x2-x21)))^2*w^2*sin(x21);
Mg3=J/r*((1-r*cos(x3))/(1+cos(x3-x22)))^2*w^2*sin(x22);
Mg4=J/r*((1-r*cos(x4))/(1+cos(x4-x23)))^2*w^2*sin(x23);
Mg5=J/r*((1-r*cos(x5))/(1+cos(x5-x24)))^2*w^2*sin(x24);
Mg6=J/r*((1-r*cos(x6))/(1+cos(x6-x25)))^2*w^2*sin(x25);
Mg7=J/r*((1-r*cos(x7))/(1+cos(x7-x26)))^2*w^2*sin(x26);
Mg8=J/r*((1-r*cos(x8))/(1+cos(x8-x27)))^2*w^2*sin(x27);
Mg9=J/r*((1-r*cos(x9))/(1+cos(x9-x28)))^2*w^2*sin(x28);
Mg10=J/r*((1-r*cos(x10))/(1+cos(x10-x29)))^2*w^2*sin(x29);
Mg11=J/r*((1-r*cos(x11))/(1+cos(x11-x30)))^2*w^2*sin(x30);
Mg12=J/r*((1-r*cos(x12))/(1+cos(x12-x31)))^2*w^2*sin(x31);
Mg13=J/r*((1-r*cos(x13))/(1+cos(x13-x32)))^2*w^2*sin(x32);
Mg14=J/r*((1-r*cos(x14))/(1+cos(x14-x33)))^2*w^2*sin(x33);
Mg15=J/r*((1-r*cos(x15))/(1+cos(x15-x34)))^2*w^2*sin(x34);
Mg16=J/r*((1-r*cos(x16))/(1+cos(x16-x35)))^2*w^2*sin(x35);
Mg17=J/r*((1-r*cos(x17))/(1+cos(x17-x36)))^2*w^2*sin(x36);
Mg18=J/r*((1-r*cos(x18))/(1+cos(x18-x37)))^2*w^2*sin(x37);
Mg19=J/r*((1-r*cos(x19))/(1+cos(x19-x38)))^2*w^2*sin(x38);
Fc1=0.5*M*Dm*((1-r*cos(x1))*(1+cos(x1-x20))^(-1))^2*w^2;
Fc2=0.5*M*Dm*((1-r*cos(x2))*(1+cos(x2-x21))^(-1))^2*w^2;
Fc3=0.5*M*Dm*((1-r*cos(x3))*(1+cos(x3-x22))^(-1))^2*w^2;
Fc4=0.5*M*Dm*((1-r*cos(x4))*(1+cos(x4-x23))^(-1))^2*w^2;
Fc5=0.5*M*Dm*((1-r*cos(x5))*(1+cos(x5-x24))^(-1))^2*w^2;
Fc6=0.5*M*Dm*((1-r*cos(x6))*(1+cos(x6-x25))^(-1))^2*w^2;
Fc7=0.5*M*Dm*((1-r*cos(x7))*(1+cos(x7-x26))^(-1))^2*w^2;
Fc8=0.5*M*Dm*((1-r*cos(x8))*(1+cos(x8-x27))^(-1))^2*w^2;
Fc9=0.5*M*Dm*((1-r*cos(x9))*(1+cos(x9-x28))^(-1))^2*w^2;
Fc10=0.5*M*Dm*((1-r*cos(x10))*(1+cos(x10-x29))^(-1))^2*w^2;
Fc11=0.5*M*Dm*((1-r*cos(x11))*(1+cos(x11-x30))^(-1))^2*w^2;
Fc12=0.5*M*Dm*((1-r*cos(x12))*(1+cos(x12-x31))^(-1))^2*w^2;
Fc13=0.5*M*Dm*((1-r*cos(x13))*(1+cos(x13-x32))^(-1))^2*w^2;
Fc14=0.5*M*Dm*((1-r*cos(x14))*(1+cos(x14-x33))^(-1))^2*w^2;
Fc15=0.5*M*Dm*((1-r*cos(x15))*(1+cos(x15-x34))^(-1))^2*w^2;
Fc16=0.5*M*Dm*((1-r*cos(x16))*(1+cos(x16-x35))^(-1))^2*w^2;
Fc17=0.5*M*Dm*((1-r*cos(x17))*(1+cos(x17-x36))^(-1))^2*w^2;
Fc18=0.5*M*Dm*((1-r*cos(x18))*(1+cos(x18-x37))^(-1))^2*w^2;
Fc19=0.5*M*Dm*((1-r*cos(x19))*(1+cos(x19-x38))^(-1))^2*w^2;
ki1=0.83;
ki2=0.83;
ki3=0.83;
ki4=0.83;
ki5=0.83;
ki6=0.83;
ki7=0.83;
ki8=0.83;
ki9=0.83;
ki10=0.83;
ki11=0.83;
ki12=0.83;
ki13=0.83;
ki14=0.83;
ki15=0.83;
ki16=0.83;
ki17=0.83;
ki18=0.83;
ki19=0.83;
ko1=0.83;
ko2=0.83;
ko3=0.83;
ko4=0.83;
ko5=0.83;
ko6=0.83;
ko7=0.83;
ko8=0.83;
ko9=0.83;
ko10=0.83;
ko11=0.83;
ko12=0.83;
ko13=0.83;
ko14=0.83;
ko15=0.83;
ko16=0.83;
ko17=0.83;
ko18=0.83;
ko19=0.83;
f(1,1)=Fa...
-(ki1*sigemai1^1.5*sin(x1)-lamda_i*Mg1/D*cos(x1))-(ki2*sigemai2^1.5*sin(x2)-lamda_i*Mg2/D*cos(x2))-(ki3*sigemai3^1.5*sin(x3)-lamda_i*Mg3/D*cos(x3))...
-(ki4*sigemai4^1.5*sin(x4)-lamda_i*Mg4/D*cos(x4))-(ki5*sigemai5^1.5*sin(x5)-lamda_i*Mg5/D*cos(x5))-(ki6*sigemai6^1.5*sin(x6)-lamda_i*Mg6/D*cos(x6))...
-(ki7*sigemai7^1.5*sin(x7)-lamda_i*Mg7/D*cos(x7))-(ki8*sigemai8^1.5*sin(x8)-lamda_i*Mg8/D*cos(x8))-(ki9*sigemai9^1.5*sin(x9)-lamda_i*Mg9/D*cos(x9))...
-(ki10*sigemai10^1.5*sin(x10)-lamda_i*Mg10/D*cos(x10))-(ki11*sigemai11^1.5*sin(x11)-lamda_i*Mg11/D*cos(x11))-(ki12*sigemai12^1.5*sin(x12)-lamda_i*Mg12/D*cos(x12))...
-(ki13*sigemai13^1.5*sin(x13)-lamda_i*Mg13/D*cos(x13))-(ki14*sigemai14^1.5*sin(x14)-lamda_i*Mg14/D*cos(x14))-(ki15*sigemai15^1.5*sin(x15)-lamda_i*Mg15/D*cos(x15))...
-(ki16*sigemai16^1.5*sin(x16)-lamda_i*Mg16/D*cos(x16))-(ki17*sigemai17^1.5*sin(x17)-lamda_i*Mg17/D*cos(x17))-(ki18*sigemai18^1.5*sin(x18)-lamda_i*Mg18/D*cos(x18))...
-(ki19*sigemai19^1.5*sin(x19)-lamda_i*Mg19/D*cos(x19));
f(2,1)=Fr...
-(ki1*sigemai1^1.5*cos(x1)+lamda_i*Mg1/D*sin(x1))*cos(fai_1)-(ki2*sigemai2^1.5*cos(x2)+lamda_i*Mg2/D*sin(x2))*cos(fai_2)...
-(ki3*sigemai3^1.5*cos(x3)+lamda_i*Mg3/D*sin(x3))*cos(fai_3)-(ki4*sigemai4^1.5*cos(x4)+lamda_i*Mg4/D*sin(x4))*cos(fai_4)...
-(ki5*sigemai5^1.5*cos(x5)+lamda_i*Mg5/D*sin(x5))*cos(fai_5)-(ki6*sigemai6^1.5*cos(x6)+lamda_i*Mg6/D*sin(x6))*cos(fai_6)...
-(ki7*sigemai7^1.5*cos(x7)+lamda_i*Mg7/D*sin(x7))*cos(fai_7)-(ki8*sigemai8^1.5*cos(x8)+lamda_i*Mg8/D*sin(x8))*cos(fai_8)...
-(ki9*sigemai9^1.5*cos(x9)+lamda_i*Mg9/D*sin(x9))*cos(fai_9)-(ki10*sigemai10^1.5*cos(x10)+lamda_i*Mg10/D*sin(x10))*cos(fai_10)...
-(ki11*sigemai11^1.5*cos(x11)+lamda_i*Mg11/D*sin(x11))*cos(fai_11)-(ki12*sigemai12^1.5*cos(x12)+lamda_i*Mg12/D*sin(x12))*cos(fai_12)...
-(ki13*sigemai13^1.5*cos(x13)+lamda_i*Mg13/D*sin(x13))*cos(fai_13)-(ki14*sigemai14^1.5*cos(x14)+lamda_i*Mg14/D*sin(x14))*cos(fai_14)...
-(ki15*sigemai15^1.5*cos(x15)+lamda_i*Mg15/D*sin(x15))*cos(fai_15)-(ki16*sigemai16^1.5*cos(x16)+lamda_i*Mg16/D*sin(x16))*cos(fai_16)...
-(ki17*sigemai17^1.5*cos(x17)+lamda_i*Mg17/D*sin(x17))*cos(fai_17)-(ki18*sigemai18^1.5*cos(x18)+lamda_i*Mg18/D*sin(x18))*cos(fai_18)...
-(ki19*sigemai19^1.5*cos(x19)+lamda_i*Mg19/D*sin(x19))*cos(fai_19);
f(3,1)=ki1*sigemai1^1.5*sin(x1)-ko1*sigemao1^1.5*sin(x20)-Mg1/D*(lamda_i*cos(x1)-lamda_o*cos(x20));
f(4,1)=ki2*sigemai2^1.5*sin(x2)-ko2*sigemao2^1.5*sin(x21)-Mg2/D*(lamda_i*cos(x2)-lamda_o*cos(x21));
f(5,1)=ki3*sigemai3^1.5*sin(x3)-ko3*sigemao3^1.5*sin(x22)-Mg3/D*(lamda_i*cos(x3)-lamda_o*cos(x22));
f(6,1)=ki4*sigemai4^1.5*sin(x4)-ko4*sigemao4^1.5*sin(x23)-Mg4/D*(lamda_i*cos(x4)-lamda_o*cos(x23));
f(7,1)=ki5*sigemai5^1.5*sin(x5)-ko5*sigemao5^1.5*sin(x24)-Mg5/D*(lamda_i*cos(x5)-lamda_o*cos(x24));
f(8,1)=ki6*sigemai6^1.5*sin(x6)-ko6*sigemao6^1.5*sin(x25)-Mg6/D*(lamda_i*cos(x6)-lamda_o*cos(x25));
f(9,1)=ki7*sigemai7^1.5*sin(x7)-ko7*sigemao7^1.5*sin(x26)-Mg7/D*(lamda_i*cos(x7)-lamda_o*cos(x26));
f(10,1)=ki8*sigemai8^1.5*sin(x8)-ko8*sigemao8^1.5*sin(x27)-Mg8/D*(lamda_i*cos(x8)-lamda_o*cos(x27));
f(11,1)=ki9*sigemai9^1.5*sin(x9)-ko9*sigemao9^1.5*sin(x28)-Mg9/D*(lamda_i*cos(x9)-lamda_o*cos(x28));
f(12,1)=ki10*sigemai10^1.5*sin(x10)-ko10*sigemao10^1.5*sin(x29)-Mg10/D*(lamda_i*cos(x10)-lamda_o*cos(x29));
f(13,1)=ki11*sigemai11^1.5*sin(x11)-ko11*sigemao11^1.5*sin(x30)-Mg11/D*(lamda_i*cos(x11)-lamda_o*cos(x30));
f(14,1)=ki12*sigemai12^1.5*sin(x12)-ko12*sigemao12^1.5*sin(x31)-Mg12/D*(lamda_i*cos(x12)-lamda_o*cos(x31));
f(15,1)=ki13*sigemai13^1.5*sin(x13)-ko13*sigemao13^1.5*sin(x32)-Mg13/D*(lamda_i*cos(x13)-lamda_o*cos(x32));
f(16,1)=ki14*sigemai14^1.5*sin(x14)-ko14*sigemao14^1.5*sin(x33)-Mg14/D*(lamda_i*cos(x14)-lamda_o*cos(x33));
f(17,1)=ki15*sigemai15^1.5*sin(x15)-ko15*sigemao15^1.5*sin(x34)-Mg15/D*(lamda_i*cos(x15)-lamda_o*cos(x34));
f(18,1)=ki16*sigemai16^1.5*sin(x16)-ko16*sigemao16^1.5*sin(x35)-Mg16/D*(lamda_i*cos(x16)-lamda_o*cos(x35));
f(19,1)=ki17*sigemai17^1.5*sin(x17)-ko17*sigemao17^1.5*sin(x36)-Mg17/D*(lamda_i*cos(x17)-lamda_o*cos(x36));
f(20,1)=ki18*sigemai18^1.5*sin(x18)-ko18*sigemao18^1.5*sin(x37)-Mg18/D*(lamda_i*cos(x18)-lamda_o*cos(x37));
f(21,1)=ki19*sigemai19^1.5*sin(x19)-ko19*sigemao19^1.5*sin(x38)-Mg19/D*(lamda_i*cos(x19)-lamda_o*cos(x38));
f(22,1)=ki1*sigemai1^1.5*cos(x1)-ko1*sigemao1^1.5*cos(x20)-Mg1/D*(lamda_i*sin(x1)-lamda_o*sin(x20))+Fc1;
f(23,1)=ki2*sigemai2^1.5*cos(x2)-ko2*sigemao2^1.5*cos(x21)-Mg2/D*(lamda_i*sin(x2)-lamda_o*sin(x21))+Fc2;
f(24,1)=ki3*sigemai3^1.5*cos(x3)-ko3*sigemao3^1.5*cos(x22)-Mg3/D*(lamda_i*sin(x3)-lamda_o*sin(x22))+Fc3;
f(25,1)=ki4*sigemai4^1.5*cos(x4)-ko4*sigemao4^1.5*cos(x23)-Mg4/D*(lamda_i*sin(x4)-lamda_o*sin(x23))+Fc4;
f(26,1)=ki5*sigemai5^1.5*cos(x5)-ko5*sigemao5^1.5*cos(x24)-Mg5/D*(lamda_i*sin(x5)-lamda_o*sin(x24))+Fc5;
f(27,1)=ki6*sigemai6^1.5*cos(x6)-ko6*sigemao6^1.5*cos(x25)-Mg6/D*(lamda_i*sin(x6)-lamda_o*sin(x25))+Fc6;
f(28,1)=ki7*sigemai7^1.5*cos(x7)-ko7*sigemao7^1.5*cos(x26)-Mg7/D*(lamda_i*sin(x7)-lamda_o*sin(x26))+Fc7;
f(29,1)=ki8*sigemai8^1.5*cos(x8)-ko8*sigemao8^1.5*cos(x27)-Mg8/D*(lamda_i*sin(x8)-lamda_o*sin(x27))+Fc8;
f(30,1)=ki9*sigemai9^1.5*cos(x9)-ko9*sigemao9^1.5*cos(x28)-Mg9/D*(lamda_i*sin(x9)-lamda_o*sin(x28))+Fc9;
f(31,1)=ki10*sigemai10^1.5*cos(x10)-ko10*sigemao10^1.5*cos(x29)-Mg10/D*(lamda_i*sin(x10)-lamda_o*sin(x29))+Fc10;
f(32,1)=ki11*sigemai11^1.5*cos(x11)-ko11*sigemao11^1.5*cos(x30)-Mg11/D*(lamda_i*sin(x11)-lamda_o*sin(x30))+Fc11;
f(33,1)=ki12*sigemai12^1.5*cos(x12)-ko12*sigemao12^1.5*cos(x31)-Mg12/D*(lamda_i*sin(x12)-lamda_o*sin(x31))+Fc12;
f(34,1)=ki13*sigemai13^1.5*cos(x13)-ko13*sigemao13^1.5*cos(x32)-Mg13/D*(lamda_i*sin(x13)-lamda_o*sin(x32))+Fc13;
f(35,1)=ki14*sigemai14^1.5*cos(x14)-ko14*sigemao14^1.5*cos(x33)-Mg14/D*(lamda_i*sin(x14)-lamda_o*sin(x33))+Fc14;
f(36,1)=ki15*sigemai15^1.5*cos(x15)-ko15*sigemao15^1.5*cos(x34)-Mg15/D*(lamda_i*sin(x15)-lamda_o*sin(x34))+Fc15;
f(37,1)=ki16*sigemai16^1.5*cos(x16)-ko16*sigemao16^1.5*cos(x35)-Mg16/D*(lamda_i*sin(x16)-lamda_o*sin(x35))+Fc16;
f(38,1)=ki17*sigemai17^1.5*cos(x17)-ko17*sigemao17^1.5*cos(x36)-Mg17/D*(lamda_i*sin(x17)-lamda_o*sin(x36))+Fc17;
f(39,1)=ki18*sigemai18^1.5*cos(x18)-ko18*sigemao18^1.5*cos(x37)-Mg18/D*(lamda_i*sin(x18)-lamda_o*sin(x37))+Fc18;
f(40,1)=ki19*sigemai19^1.5*cos(x19)-ko19*sigemao19^1.5*cos(x38)-Mg19/D*(lamda_i*sin(x19)-lamda_o*sin(x38))+Fc19;
var1=[0.53 0.53 0.53 0.53 0.53 0.53 0.53 0.53 0.53 0.53 0.53 0.53 0.53 0.53 0.53 0.53 0.53 0.53 0.53 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0.37 0 0.0005]';
tol=0.001;
nmax=500; 
numvar=length(var); 
numeq=length(f); 
for j=1:numeq 
for i=1:numvar
    J(j,i)=diff(f(j),var(i)); 
end 
end 
n=0;  
dmain=1.1*tol*ones(size(var,1)); 
while(any(abs(dmain)>tol)&&(n<nmax)) 
    f1=subs(f,var,var1); 
    J1=subs(J,var,var1); 
    if(abs(det(J1))>tol*tol) 
    dmain=inv(J1)*f1; 
    end
    var1=var1-dmain; 
    n=n+1; 
end 
if(n==nmax) 
    error('no feasible solution'); 
end 
return

Best Answer

Your line

J=0.4*M*(D/2)^2;

implies that J must be numeric. But later you have

for j=1:numeq 
for i=1:numvar
    J(j,i)=diff(f(j),var(i)); 
end 
end

which is a symbolic computation being assigned into that numeric variable.

Your var1 needs to be a row vector not a column vector at the time of the subs(), so

    f1 = subs(f,var,var1.'); 
    J1 = subs(J,var,var1.');

You should probably be taking double(f1) and double(J1) for your computations.

Instead of using inv(J1)*f1 you should be using J1\f1

Putting these together,

while(any(abs(dmain)>tol)&&(n<nmax)) 
  f1 = double( subs(f,var,var1.') ); 
  J1 = double( subs(J,var,var1.') ); 
  if (abs(det(J1))>tol*tol) 
    dmain = J1 \ f1; 
  end
  var1 = var1 - dmain; 
  n=n+1; 
end

Related Solutions

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

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.

MATLAB: Fsolve doesn’t give correct answer even with close initial guess

A serious issue is that you use many numbers accurate to only 2 decimal digits. So even though you may think the answer is correct, it can easily be not a solutrion at all. dpb showed that to be a fact. As well, any time you have multiple trig functions in there, expect any solution to be non-unique.

Let me try your start points. I'll call x0 the value that you CLAIM to be a solution. x0 is the vector that you started fsolve at.

x0=[2.55707, 3.13087,0.68765,0,2.443221,2.557073];
>> xinit=[2.5570, 3.005,1,0,2.3685,2.0853];
>> [IK(xinit),IK(x0)]
ans =
      0.86398            1
      0.47215       -3e-06
     -0.17497  -8.3398e-07
     -0.40696        3e-06
      0.85941            1
      0.30952    5.518e-07
       0.2965   8.3398e-07
     -0.19621   -5.518e-07
      0.93466            1
      0.97211      0.81522
     -0.77106     -0.66715
      0.68671       1.0763
>> norm(IK(x0))
ans =
       2.2952
>> norm(IK(xinit))
ans =
       2.2385

As you can see, evaluating IK at each of those points gives very different results, yet the norm of the objective is quite close. That bodes very poorly for your claim that this is truly a valid solution.

What does fsolve give us?

[xf,fv,exitflag] = fsolve(@IK,x0)
Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using Levenberg-Marquardt algorithm instead. 
> In fsolve (line 310) 
No solution found.
fsolve stopped because the last step was ineffective. However, the vector of function
values is not near zero, as measured by the value of the function tolerance. 
<stopping criteria details>
xf =
       4.1679      0.52101     -0.77606      0.49216       1.4093       4.2109
fv =
      0.93396
      0.12102
      0.33626
   -0.0022887
      0.94293
     -0.33299
     -0.35737
      0.31023
      0.88093
   5.0244e-06
  -2.8169e-06
   1.2031e-08
exitflag =
    -2

In fact, fsolve is not terribly happy with the result. HOWEVER, it is considerably better than either your start point OR your claimed true "solution".

norm(IK(xf))
ans =
       1.7321

So fsolve did indeed improve the result. With 12 equations in 6 unknowns, expect there to be no true solution. Worse, when you use 2 digit numbers in such a complex problem, expect garbage out the end. The set in xf is indeed an improvement, by any measure we can apply.

[IK(xinit),IK(x0),IK(xf)]
ans =
      0.86398            1      0.93396
      0.47215       -3e-06      0.12102
     -0.17497  -8.3398e-07      0.33626
     -0.40696        3e-06   -0.0022887
      0.85941            1      0.94293
      0.30952    5.518e-07     -0.33299
       0.2965   8.3398e-07     -0.35737
     -0.19621   -5.518e-07      0.31023
      0.93466            1      0.88093
      0.97211      0.81522   5.0244e-06
     -0.77106     -0.66715  -2.8169e-06
      0.68671       1.0763   1.2031e-08

If you want a better solution, then I would start by using more accurate coefficients inside the objective. Then I would recognize that the solution will not be unique, and that with more equations than unknowns, the solution will generally not be exact either.

Best Answer

Related Solutions

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

MATLAB: Fsolve doesn’t give correct answer even with close initial guess

Related Question