论坛 › Matlab讨论区 › fsolve求解非线性方程组时遇到的问题

buaa_32051114 发表于 2007-8-10 12:43

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

我是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到u是挠度，应该对称。可以用来检验结果是否正确。

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

程序：
functionwork

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)
   

    functionf = 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);
                  elseif80e9*(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 发表于 2007-8-10 12:45

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

buaa_32051114 发表于 2007-8-12 21:20

自己解决了，呵呵。

无水1324 发表于 2007-8-13 09:31

回复 #3 buaa_32051114 的帖子

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

hoozecn 发表于 2007-8-13 15:20

都是if else 啊 呵呵～～

查看完整版本: fsolve求解非线性方程组时遇到的问题