最近在写论文，需要求解一个由十二个方程组成的超越方程组，试用newton-Raphon公式求解过，但对初值的选取一直很烦恼，还望各位高手能不吝赐教，万分感谢！
function f=myfun(x)
c(1)=x(1);
c(2)=x(2);
k(1)=x(3);
k(2)=x(4);
k(3)=x(5);
k(4)=x(6);
k(5)=x(7);
k(6)=x(8);
k(7)=x(9);
k(8)=x(10);
A=x(11);
B=x(12);
%A处位移连续
f(1)=exp(A*0.9931)*(c(1)*cos(A*0.9931)+c(2)*sin(A*0.9931))-(-A^4*8751.6/24+k(1)*A^3/6+k(2)*A^2/2+k(3)*A+k(4))/(1.1248*10^6);
%A处位移的一阶导数连续
f(2)=exp(A*0.9931)*0.9931*((c(1)*(cos(A*0.9931)-sin(A*0.9931))+c(2)*(cos(A*0.9931)+sin(A*0.9931)))-(-A^3*8751.6/6+k(1)*A^2/2+k(2)*A+k(3))/(1.1248*10^6));
%A处位移的二阶导数连续
f(3)=2*exp(A*0.9931)*0.9931^2*(-c(1)*sin(A*0.9931)+c(2)*cos(A*0.9931))-(-A^2*8751.6/2+k(1)*A+k(2))/(1.1248*10^6);
%A处位移的三阶导数连续
f(4)=2*exp(A*0.9931)*0.9931^3*((c(1)*(-cos(A*0.9931)-sin(A*0.9931))+c(2)*(cos(A*0.9931)-sin(A*0.9931)))-(-A*8751.6+k(1))/(1.1248*10^6));
%A处位移的为D(s)
f(5)=(-A^4*8751.6/24+k(1)*A^3/6+k(2)*A^2/2+k(3)*A+k(4))/(1.1248*10^6)-0.008;
%B处位移连续
f(6)=(B^4*8751.6/24+k(1)*B^3/6+k(2)*B^2/2+k(3)*B+k(4))/(1.1248*10^6)-exp(B*0.9931)*(k(5)*cos(B*0.9931)+k(6)*sin(B*0.9931))-exp(-B*0.9931)*(k(7)*cos(B*0.9931)+k(8)*sin(B*0.9931))-10;            
%B处位移的一阶导数连续
f(7)=(B^3*8751.6/6+k(1)*B^2/2+k(2)*B+k(3))/(1.1248*10^6)-exp(B*0.9931)*0.9931*(k(5)*(cos(B*0.9931)-sin(B*0.9931))+k(6)*(cos(B*0.9931)+sin(B*0.9931)))+exp(-B*0.9931)*0.9931*(k(7)*(cos(B*0.9931)+sin(B*0.9931))+k(8)*(-cos(B*0.9931)+sin(B*0.9931)));
%B处位移的二阶导数连续
f(8)=(B^2*8751.6/2+k(1)*B+k(2))/(1.1248*10^6)-2*exp(B*0.9931)*0.9931^2*(-k(5)*sin(B*0.9931)+k(6)*cos(B*0.9931))-2*exp(-B*0.9931)*0.9931^2*(k(7)*sin(B*0.9931)-k(8)*cos(B*0.9931));
%B处位移的三阶导数连续
f(9)=(B*8751.6+k(1))/(1.1248*10^6)-2*exp(B*0.9931)*0.9931^3*(k(5)*(-cos(B*0.9931)-sin(B*0.9931))+k(6)*(cos(B*0.9931)-sin(B*0.9931)))-2*exp(-B*0.9931)*0.9931^3*(k(7)*(cos(B*0.9931)-sin(B*0.9931))+k(8)*(cos(B*0.9931)+sin(B*0.9931)));
%B处位移的为(delta-D(s))
f(10)=(B^4*8751.6/24+k(1)*B^3/6+k(2)*B^2/2+k(3)*B+k(4))/(1.1248*10^6)-(10-0.008);
%滑坡中部W/2处转角为零
f(11)=exp(40*0.9931)*(k(5)*(cos(40*0.9931)-sin(40*0.9931))+k(6)*(cos(40*0.9931)+sin(40*0.9931)))+exp(-40*0.9931)*(k(7)*(-cos(40*0.9931)-sin(40*0.9931))+k(8)*(cos(40*0.9931)-sin(40*0.9931)));
%整体力平衡条件
f(12)=exp(A*0.9931)*(c(1)*(cos(A*0.9931)+sin(A*0.9931))+c(2)*(-cos(A*0.9931)+sin(A*0.9931)))...
    +exp(B*0.9931)*(k(5)*(-cos(B*0.9931)-sin(B*0.9931))+k(6)*(cos(B*0.9931)-sin(B*0.9931)))...
    +exp(-B*0.9931)*(k(7)*(cos(B*0.9931)-sin(B*0.9931))-k(8)*(-cos(B*0.9931)-sin(B*0.9931)))...
    +exp(40*0.9931)*(k(5)*(cos(40*0.9931)+sin(40*0.9931))+k(6)*(-cos(B*0.9931)+sin(B*0.9931)))...
    +exp(-40*0.9931)*(k(7)*(-cos(40*0.9931)+sin(40*0.9931))-k(8)*(-cos(B*0.9931)-sin(B*0.9931)))...
    -(0.008*(A+B)*2*0.9931);
f=[f(1) f(2) f(3) f(4) f(5) f(6) f(7) f(8) f(9) f(10) f(11) f(12)]

[[it] 本帖最后由 piaoyun003 于 2008-9-24 11:17 编辑 [/it]]