作者trashken (上个昵称)
看板Fortran
标题[请益] 与积分有关系得非线性方程
时间Thu Sep 8 16:52:34 2011
程式原始码:
PROGRAM guess
use IMSL
IMPLICIT none
external FCN
real, parameter :: ERRREL=0.0001
integer, parameter :: N=1
integer, parameter :: ITMAX=100
real :: XGUESS(N)=(/0.4/)
real X(N), FNORM
CALL NEQNF(FCN, ERRREL, N, ITMAX, XGUESS, X, FNORM)
write(*,*) X
stop
end
subroutine FCN (XA, F, N)
implicit none
integer N
real, target ::XA(N)
real F(N)
real*8 ::El,kT,ul,ur,Tl,Tr,pi,fl,fr,Ea,T,VL,Y
REAL,POINTER ::N0
N0=>XA(1)
El=450
kT=0
ul=60
ur=60
VL=60
Tl=1
Tr=1
T=Tl+Tr
pi=3.14159
fl=1/(exp((Ea-ul)/kT)+1)
fr=1/(exp((Ea-ur)/kT)+1)
do Ea=0.d0,300.d0,0.1d0
Y=Y+(-(fl*Tl+fr*Tr))/(pi*(Tl+Tr))*((-(1-N0)*(T/2))/((Ea-El)**2+(T/2)**2)-(N0*T
/2)/((Ea-El-VL)**2+(T/2)**2))
end do
F(1)=Y-N0
RETURN
END Subroutine
=============================================================================
可不可以请问一下版上各位大大
我目前在解一个与积分有关的方程式
也可以说是一个耦合方程式
要积分的式子是
(-(fl*Tl+fr*Tr))/(pi*(Tl+Tr))*((-(1-N0)*(T/2))/((Ea-El)**2+(T/2)**2)-(N0*T
/2)/((Ea-El-VL)**2+(T/2)**2))
而积分的变数是Ea 而要解的未知数是N0
而N0也等於这一长串要积分的式子
也就是说N0=(-(fl*Tl+fr*Tr))/(pi*(Tl+Tr))*((-(1-N0)*(T/2))/((Ea-El)**2+(T/2)**2
)-(N0*T/2)/((Ea-El-VL)**2+(T/2)**2))
而我使用函式库的NEQNF功能
而使用一个简单的DO回圈来积那一长串的式子
再令这一长串的式子为Y
而後在DO回圈外使用"F(1)=Y-N0"这个式子 来求解这个未知数N0
程式可以跑
但问题是跑出来的值跟我要模拟的图相比较之下都怪怪的
在这里可否请教各位 如果欲求解的未知数出现在积分式中
那要怎麽样来写方程式 才会比较恰当
我这样用DO回圈直接求解的方式 好像错了
这个问题已经卡住小弟两三个月了
希望有大大可以替小弟解惑
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 140.115.72.52
※ 编辑: trashken 来自: 140.115.72.52 (09/08 16:54)
1F：推 terryys:比较好的做法是把函数另外写作一个副程式,用imsl之类的 09/08 17:44
2F：→ terryys:算积分 你的FCN至少有三个错 09/08 17:45
3F：推 terryys:fl fr随Ea变动,应该把他们放到do里面 09/08 17:49
4F：→ terryys:do Ea=0.d0,300.d0,0.1d0 这个写法已经过时,应该用整数 09/08 17:50
5F：→ trashken:可以请教一下我FCN的错在哪里吗 还有IMSL计算积分的话 09/08 17:52
6F：→ trashken:大概要用什麽函式库 09/08 17:53
7F：→ terryys:然後算积分的时候应该还要乘上dEa,你这里就是0.1 09/08 17:53
8F：→ trashken:是整个式子乘上0.1吗 09/08 17:54
9F：→ trashken:我手上有彭国伦写的fortran参考书 可以请教一下要写积分 09/08 18:04
10F：→ trashken:副程式 怎麽样比较快 09/08 18:04
11F：→ trashken:我乘上0.1了 结果还是没变 09/08 18:12
12F：推 terryys:我上面说的那两个改了吗? 还有Y没有初始化,虽然应该都是0 09/08 20:13
13F：→ terryys:但是有初始化比较心安理得 09/08 20:13
14F：推 terryys:imsl我没有用过，但是google很容易就可以找到应该用哪一个 09/08 20:24
15F：推 s06yji3:我之前有解过类似的题目，我用Simpson3-8搭配NR 09/10 17:14
16F：→ s06yji3:Simpon3-8和NR我是另外写个副程式丢进去算 09/10 17:15