作者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