fsolve求解非线性方程组时遇到的问题

buaa_32051114

我是matlab的新手，前几天要编个程序需要用matlab先算出初值。不过结果始终不满意，我先简单介绍下问题，再给出自己程序，希望高手们帮我看看。谢谢了！
  
          位移u(数组),然后是中间量a,b,c （也都是数组）是u 的函数。根据这些算出应力q1,q2,q3..再是根据应力算出力和力偶n,m并得到加速度。最后是用a,b,c,m,n得到下一时刻位移u.  
近似后，u = (1/4 )*a*t^2,作为方程。

这是梁的位移，从u[21]到u[41]是挠度，应该对称。可以用来检验结果是否正确。

我是新手，可能是最基本的地方错了，具体的某个量应该没什么问题，如果不理解可以不用考虑。

程序：
function  work

format long
u0 = zeros(42,1);
u = zeros(42,1);
a = zeros(21,1);
b = zeros(21,1) ;
c = zeros(21,1);
q1 = zeros(21,5);
p2 = zeros(21,5);
q3 = zeros(21,1);
n = zeros(21,1);
m = zeros(21,1);



u = fsolve(@myfun,u0)
   

    function  f = myfun(u)
        u(1) = 0; u(21) = 0; u(22) = 0; u(42) = 0;

        
            a(1) = u(2)/0.01;
            a(21) = -u(20)/0.01;
            for i=2:20
                a(i)=(u(i+1)-u(i-1))/(2*0.01);
            end
        
     
            b(1) = u(23)/0.01;
            b(21) = -u(41)/0.01;
            for i=2:20
                b(i)=(u(i+22)-u(i+20))/(2*0.01);
            end

        
         
            c(1) = (u(24)-2*u(23))/(0.01)^2;
            c(21) = (u(40)-2*u(41))/(0.01)^2;
            for i=2:20
                c(i) = (u(i+22)-2*u(i+21)+u(i+20))/(0.01)^2;
            end
   
        
        
            for i=1:21
                for j=1:5
                    if abs(80e9*(a(i)-j*(4e-4)*c(i)+0.5*b(i)^2))<=0.3e9
                        q1(i,j)=80e9*(a(i)-j*(4e-4)*c(i)+0.5*b(i)^2);
                    elseif 80e9*(a(i)-j*(4e-4)*c(i)+0.5*b(i)^2)>0.3e9
                            q1(i,j)=0.3e9;
                    else
                            q1(i,j)=-0.3e9;
                    end
                end
            end

      
         
            for i=1:21
                for j=1:5
                    if abs(80e9*(a(i)+j*(4e-4)*c(i)+0.5*b(i)^2))<=0.3e9
                        p2(i,j)=80e9*(a(i)+j*(4e-4)*c(i)+0.5*b(i)^2);
                    elseif  80e9*(a(i)+j*(4e-4)*c(i)+0.5*b(i)^2)>0.3e9
                            p2(i,j)=0.3e9;
                    else
                            p2(i,j)=-0.3e9;
                    end
                end
            end

            
         
            for i=1:21
               if abs(80e9*(a(i)+0.5*b(i)^2))<=0.3e9
                    q3(i)=80e9*(a(i)+0.5*b(i)^2);
                elseif 80e9*(a(i)+0.5*b(i)^2)>0.3e9
                    q3(i)=0.3e9;
                else
                    q3(i)=-0.3e9;
                end
            end
            
            
         
        
            for i=1:21
                n(i) = (q1(i,5) + p2(i,5))*0.02*(4e-4)/2+q3(i)*0.02*(4e-4);
                for j=1:4
                    n(i) = n(i) + (q1(i,j)+p2(i,j))*0.02*(4e-4);
                end
            end
        
        
        
            for i = 1:21
                m(i) = -(q1(i,5)-p2(i,5))*0.02*5*(4e-4)^2/2;
                for j=1:4
                    m(i) = m(i) -(q1(i,j)-p2(i,j))*j*0.02*(4e-4)^2;
                end
            end

  
  for i=2:20
         f(i) = (1/4)*(((1+a(i+1))*n(i+1) - (1+a(i-1))*n(i-1) + m(i+1)*c(i+1) - m(i-1)*c(i-1))/(0.01*2*0.02*(4e-3)*2700)         +b(i)*96500/((4e-3)*2700))*(5e-4)^2 -u(i);
     end
     for i=23:41
         f(i) = (1/4)*(-(m(i-20) - 2*m(i-21) + m(i-22)/(0.02*(4e-3)*2700*(0.01)^2)) + 96500/((4e-3)*2700))*(5e-4)^2 - u(i);
     end

buaa_32051114

上次那个发的时候出了点问题，重发一个。不好意思:@L

buaa_32051114

自己解决了，呵呵。

无水1324

能否介绍一下解决的方法及最后的结论，供后来者参考，谢谢！

hoozecn

都是if else 啊 呵呵～～