| 查看: 814 | 回复: 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 |
» 猜你喜欢
拟解决的关键科学问题还要不要写
已经有6人回复
存款400万可以在学校里躺平吗
已经有14人回复
Materials Today Chemistry审稿周期
已经有6人回复
基金委咋了?2026年的指南还没有出来?
已经有10人回复
基金申报
已经有6人回复
推荐一本书
已经有13人回复
国自然申请面上模板最新2026版出了吗?
已经有17人回复
纳米粒子粒径的测量
已经有8人回复
疑惑?
已经有5人回复
计算机、0854电子信息(085401-058412)调剂
已经有5人回复
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











回复此楼