声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1338|回复: 4

[综合讨论] fsolve求解非线性方程组时遇到的问题

[复制链接]
发表于 2007-8-10 12:43 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
我是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
回复
分享到:

使用道具 举报

 楼主| 发表于 2007-8-10 12:45 | 显示全部楼层
上次那个发的时候出了点问题,重发一个。不好意思:@L
 楼主| 发表于 2007-8-12 21:20 | 显示全部楼层
自己解决了,呵呵。
发表于 2007-8-13 09:31 | 显示全部楼层

回复 #3 buaa_32051114 的帖子

能否介绍一下解决的方法及最后的结论,供后来者参考,谢谢!
发表于 2007-8-13 15:20 | 显示全部楼层
都是if else 啊 呵呵~~
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-24 21:41 , Processed in 0.054229 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表