| 查看: 869 | 回复: 2 | ||
[求助]
fortran程序总是出现“root must be bracketed for zbrent”求指点
|
|
-13.34=2logA+3logB+2(-6.6B+0.045A-0.001A2)+3(-0.2B-3.9A+1.7A2) 这个是我想解的方程，想输入一个A就能对应得到一个B A取值范围是10e-5-5e-2，B取值范围是5e-6--1e-4 希望高手能帮忙解答一下，小弟感激不尽，下面是程序 PROGRAM main INCLUDE 'DEFIN.DIM' EXTERNAL func x1=1.D-5 x2=5.D-2 tol=1.D-10 999 BB=zbrent(x1,x2,tol) GOTO 999 WRITE(*,*)BB END FUNCTION zbrent(x1,x2,tol) INCLUDE 'DEFIN.DIM' c INTEGER ITMAX c REAL zbrent,tol,x1,x2,func,EPS EXTERNAL func PARAMETER (ITMAX=100,EPS=3.e-8) c INTEGER iter c REAL a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm a=x1 b=x2 fa=func(a) fb=func(b) if((fa.gt.0..and.fb.gt.0.).or.(fa.lt.0..and.fb.lt.0.))pause *'root must be bracketed for zbrent' c=b fc=fb do 11 iter=1,ITMAX if((fb.gt.0..and.fc.gt.0.).or.(fb.lt.0..and.fc.lt.0.))then c=a fc=fa d=b-a e=d endif if(abs(fc).lt.abs(fb)) then a=b b=c c=a fa=fb fb=fc fc=fa endif tol1=2.*EPS*abs(b)+0.5*tol xm=.5*(c-b) if(abs(xm).le.tol1 .or. fb.eq.0.)then zbrent=b goto 99 c return endif if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then s=fb/fa if(a.eq.c) then p=2.*xm*s q=1.-s else q=fa/fc r=fb/fc p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.)) q=(q-1.)*(r-1.)*(s-1.) endif if(p.gt.0.) q=-q p=abs(p) if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then e=d d=p/q else d=xm e=d endif else d=xm e=d endif a=b fa=fb if(abs(d) .gt. tol1) then b=b+d else b=b+sign(tol1,xm) endif fb=func(b) 11 continue pause 'zbrent exceeding maximum iterations' zbrent=b 99 return end FUNCTION func(BB) INCLUDE 'DEFIN.DIM' AA=1.D-5 func=2*log(AA)+3*log(BB)+2*(-6.6*BB+0.045*AA-0.001*AA**2) + +3*(-0.2*BB-3.9*AA+1.7*AA**2)+13.34 c END IF RETURN END |
回复此楼» 猜你喜欢
河北省自然科学基金
已经有8人回复
-
西安交大新媒学院副院长用撤稿论文结题
已经有5人回复
-
论文撤稿了
已经有5人回复
-
某211大学教师把个人教师官方主页改成:我跑了我跑了我跑了!官宣跑路！
已经有5人回复
-
26/27申博自荐
已经有9人回复
-
售SCI一区T0P文章，我:8.O.5.5.1.O.5.4,科目齐全,可+急
已经有3人回复
-
售SCI一区T0P文章，我:8.O.5.5.1.O.5.4,科目齐全,可+急
已经有7人回复
-
揭秘青基评审内幕：几个A才能顺利中标
已经有4人回复
-
青B发送上会通知了吗
已经有7人回复
-
博士申请
已经有3人回复
baobiao007
木虫 (职业作家)
中国特色
- 应助: 201 (大学生)
- 金币: 6482.7
- 散金: 557
- 红花: 40
- 帖子: 3050
- 在线: 1009.9小时
- 虫号: 505962
- 注册: 2008-02-18
- 专业: 应用地球物理学
2楼2013-10-25 10:52:23
3楼2013-10-25 13:14:30
10