[請益] 與積分有關係得非線性方程

[trashken (上個暱稱)]

程式原始碼:

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)
推 terryys:比較好的做法是把函數另外寫作一個副程式,用imsl之類的
→ terryys:算積分  你的FCN至少有三個錯
推 terryys:fl fr隨Ea變動,應該把他們放到do裡面
→ terryys:do Ea=0.d0,300.d0,0.1d0  這個寫法已經過時,應該用整數
→ trashken:可以請教一下我FCN的錯在哪裡嗎 還有IMSL計算積分的話
→ trashken:大概要用什麼函式庫
→ terryys:然後算積分的時候應該還要乘上dEa,你這裡就是0.1
→ trashken:是整個式子乘上0.1嗎
→ trashken:我手上有彭國倫寫的fortran參考書 可以請教一下要寫積分
→ trashken:副程式 怎麼樣比較快
→ trashken:我乘上0.1了 結果還是沒變
推 terryys:我上面說的那兩個改了嗎? 還有Y沒有初始化,雖然應該都是0
→ terryys:但是有初始化比較心安理得
推 terryys:imsl我沒有用過，但是google很容易就可以找到應該用哪一個
推 s06yji3:我之前有解過類似的題目，我用Simpson3-8搭配NR
→ s06yji3:Simpon3-8和NR我是另外寫個副程式丟進去算