ÂùÂIÃä¬É­È°ÝÃD¡G®gÀ»ªk

 

®gÀ»ªk

§Ú­Ì«e­±¤w¸g´£¨ì®gÀ»ªkªº°ò¥»µ¦²¤¡A´N¬O¦b°_ÂIªº¨º¤@°¼°w¹ï¤£¨¬ªº±ø¥ó¥[¥H²q´ú¡A¦p¦¹±o¥H§Q¥Î¼Ð·Çªº¡]­nµ¹ªì©l­Èªº¡^ºtºâªk¦p Runge-Kutta ¨Ó§â¿n¤À¤@¸ô°µ¨ì¥t¤@°¼²×ÂI¡C ³o¼Ë©Ò±o¥X¨Óªº¸Ñ·íµM¤£·|«ê¦n§k¦X¥t¤@°¼ªºÃä¬É±ø¥ó¡A¦ý¬O§Ú­Ì§Æ±æ§Q¥Î¦¹¤@½üªºµ²ªG¡A¨Ó¬°¤U¤@¦¸·sªº²q´ú´£¨Ñ§ó¦nªº«Øij¡Aª½¨ì³Ì«á²×¯à©R¤¤¥Ø¼Ð¡]¤]´N¬Oº¡¨¬¤FÃä¬É±ø¥ó¦Ó§¹¦¨¤F°ÝÃDªº¨D¸Ñ¡C

¦]¦¹¡A¨D¸Ñ°ÝÃDªºÃöÁä¦b©ó¡A¦p¦ó«Ø¥ß©Ò»Ýªº¾÷¨î¨ÓÅý§Ú­Ìªºªì©l²q´úªº©R¤¤«×¡A±o¥H³z¹L­¡¥Nªº¤è¦¡¤£Â_§ï¶i¡H

¤@­ÓÁo©úªº¿ìªk¬O§â³o­Ó§Æ±æ¤£Â_´£ª@©R¤¤«×ªº°Ê§@Âà¤Æ¬°¬O¨D®Úªº°ÝÃD¡C§Ú­Ì¥i¥H¥O F(V) ³o¼Ë¤@­Ó¥H V ¬°¦ÛÅܼƪº¨ç¼Æ¡A¥¦¬O¦b²q´úªì©l±ø¥ó V ©Ò°µ¥X¨Óªº y(x2;V) »PÀ³¸Óº¡¨¬ªºÃä¬É±ø¥ó¤§¶¡ªº®t¡A§Y F(V) = y(x2;V) - y(x2;BC)¡CÅã¦Ó©ö¨£¦a¡A­Y¤w©R¤¤Ãä¬É±ø¥ó¡A«h F(V) = 0¡A§_«h F(V) ´N¤£µ¥©ó¹s¡C¤@¥¹«Ø¥ß¤F F(V) ³o¼Ëªº¨ç¼Æ¡A§Ú­Ì´N¥i¥H¥Î¼Ð·Çªº¤èªk¨D®Ú¡A¨D¨ì¤F F(V) ªº®Ú¡A¾ã­ÓÂùÂIÃä¬É­Èªº°ÝÃD¤]´N¨D¸Ñ§¹¦¨¡C

¦b¦¹¡A§Ú­Ì»Ý­n¨D¸Ñ V¡A§Q¥Î¤û¹yªk¡A¬O¬Û·í©ó¤@¨B¨B¨«¦V

Vnew = Vold + dV

¨ä¤¤ dV º¡¨¬¤U­±Ãö«Y¦¡

J¡DdV = -F

³oùئ³¯A¤Î¨ì«e­±¥X²{ªº Jacobian Matrix

«ÜÅãµM¦a¡A§Ú­Ì®Ú¥»¤£ª¾¹D F(V) ªº¨ãÅé¨ç¼Æ«¬¦¡¡A¦]¦¹µL±q°µ¸ÑªR°¾·L¤À¡A¥u¯à«ö·L¤Àªº°ò¥»©w¸q¡A¥Ñ¥H¤U¼Æ­È·L¤Àªº¤è¦¡Àò±o±o Jij¡G

¦³¤F J ¤§«á¡A¦ÛµM¥i¥Î¸Ñ½u©Ê¥N¼Æ¤èµ{¦¡ªº¤èªk¨Ó¨D¥X dV ¡C

³oùØÁÙ¦³¤@­ÓÃB¥~ªºÂI­n¦Ò¶q¡A¯Âºéªº¤û¹yªk¥»¨­¦³¤@¨Ç¯ÊÂI¡A§Q¥Î©Ò¿× "¥þ°ì¦¬ÀĪº¤û¹yªk"¡A¥i¤j¤j´£ª@§ä¨ì®Úªº¯à¤O¡A°ò¥»¤W¥¦¬O¤û¹yªk¦A»²¥Hºë¥©ªº·j´M¾÷¨î¡A¸Ô±¡½Ð¨£¬ÛÃö³¹¸`¡C

 

 

 

°Æµ{¦¡ªº¨Ï¥Î

¥»³æ¤¸©Ò·|¥Î¨ìªº°Æµ{¦¡¡A¤£§t¨Ï¥ÎªÌ©Ò»Ý­n¼gªº¡A¤@¦@¦³¤Q­Ó¤§¦h¡A½Ð¤j®aª`·N¡C¡]©Ò¦³ newt ¬ÛÃö»P odeint ¬ÛÃöªº³£­n¥Î¨ì¡AùØÀY¥]§t¡G¥þ°ì¦¬ÀĤû¹yªk»Ý¥Î newt¡Blnsrch¡Bfmin¡Bfdjac¡Bludcmp¡Blubksb ¡F¨B´T½Õ±±¶©¥¨®w¶ðªk»Ý¥Î odeint¡B rkqs¡Brkck¡F®gÀ»ªk»Ý¥Î shoot¡C¡^

¥»¸`ªº¥D­n°Æµ{¦¡¬O shoot ¡A¥¦·|¶Ç¤Jªì©l²q´ú­È¡A¨ÃÅX°Ê rkqs ¡]¦³½Õ¾A¨B´Tªº Runge-Kutta ¸Ñ±`·L¤è¤èªk¡^¨«§¹¤@±ø¸Ñ¡C«Ü­«­nªº¬O¥¦¥²¶·³Q¡]¥þ°ì¦¬ÀĪº¡^¤û¹yªk newt ¥s¥Î¡]§Ú­Ì«e­±¦³¤¶²Ð¹L¡Anewt ªº¥Îªk¬Û·í³æ¯Â¡^¡A¦]¦¹½Ðª`·N¨ì³oùØ shoot.f Àɤº¹ê»Úªº°Æµ{¦¡¦W¬O funcv ¡C

(°Æµ{¦¡)
¨Ï¥ÎªÌ¥²¶·¦Û¦æ´£¨Ñªº°Æµ{¦¡¦b³oùØ«h¦³¤T­Ó¡A¤À§O¬O (1) load¡B(2) score¡B(3) derivs¡C¦A¦¸´£¿ô¡Ashoot.f °Æµ{¦¡¬O¥s§@ funcv¡A³o¦¸ªº funcv ¤£¥²¨Ï¥ÎªÌ¦Û¤v¼g¡C

(¤½¦@°Ï¶ô)
¼¶¼g¥Dµ{¦¡®É¡A­nª`·N common block ªº³]©w¡A¥Ñ shoot ªº¤º®e¥i¥Hª¾¹D¡A¦³¨ÇÅܼƬO­n³z¹L common blcok ¶Çµ¹°Æµ{¦¡ªº¡A¥t¥~¡A¦pªG­nµe¥X¦±½u¡A«h rkqs ªº¤¤¶¡¹Lµ{µ²ªG¤]­n¥Ñ common blcok ¨ú¥X¨Ó¡C½Ðª`·N½Ò¥»ªþªº shoot (funcv) °Æµ{¦¡·|§â kmax ³]¦¨ 0 ¡A¦]¦¹¤£·|§@¤¤¶¡¹Lµ{ªºÀx¦s¡A·Q¥Î¨ì³o¨Ç¤¤¶¡µ²ªG´N­n¦Û¦æ§â shoot ¤¤¤§ kmax §ï³]¡]³Ì¦h¤£¯à¶W¹L KMAXX¡A§@ªk°Ñ¦Ò odeint ¤§»¡©ú¡^¡A©Î¬O§Q¥Î¨Ï¥ÎªÌ¦Û¤v¼gªº°Æµ{¦¡¤¤¦s¨ú common block ¡A¦b©I¥s odeint ¤§«e´N¥ý§â kmax ¦A§ï¦¨«D¹sªº­È¡C

(parameter »P common)
¥t¥~¡A«Å§i¤F parameter ªº§ô¦è¡A¤]¤£¯à©ñ¦b common ¤¤¡A¤ñ¤è»¡¡Anvar ¦³³Q«Å§i¦b common block ¤¤¡Anvar ´N¤£¯à³Q©w¬° parameter¡A¦³»Ý­n§@ nvar ¤j¤p°}¨ì«Å§iªº¦a¤è¡A´N­n¥t¦æ¥Î¨ä¥L°Ñ¼Æ©Îª½±µµ¹¼Æ¦r¨Ó«Å§i¡C ­È±o¤@´£ªº¬O¡Aodeint ¤º³¡Áö¦³«Å§i¥Î¨ì ystart(nvar) ³o¼ËªºªF¦è¡A¦ý ystart ¬O¥Ñ¥~³¡ y ¶Ç¶i¥hµ¹ odeint ªº¡A¥B nvar ¤]¦Û¥~³¡¶Ç¤J¡A«ö Fortran ¦V°Æµ{¦¡¶Ç»¼°}¦Cªº¤è¦¡¡A¬O¥u¶Ç°}¦Cªº°_©l¦ì¸m¡A¦]¦¹¦b odeint °Æµ{¦¡¤¤¨Ã¨S¦³·s¶}°O¾ÐÅé¶J¦sªºªÅ¶¡¡A§Ú­Ì¦]¦¹¤£¥²§â nvar ¦b¤W¼hµ{¦¡³]©w¦¨±`¼Æ¡A²¦³º¤£·|¯uªº¦b odeint ¤¤¦³°t¸m nvar ¤j¤p¤§°}¦C³oºØ°Ê§@¡A¨Æ¹ê¤W¡A¦¹¨t¦C°Æµ{¦¡®M²Õ¤v¹w³] nvar ªº¤W­­¬° NMAX¡A¤]´N¬O 50¡A¤wµ´¹ï°÷¥Î¤F¡C

(¿n¤Àªì©l­È)
¨Ï¥ÎªÌ©Ò­n¼gªº¥Dµ{¦¡¡A¥¦·|©I¥s newt¡A ¥Dµ{¦¡­n¥ýŪ¤Jªì©l²q´úªº V(n2)¡A¶Çµ¹ newt ·í§@¬O­n¨Dªº®Ú¡AµM«á newt ¶Çµ¹ shoot (funcv)¡A¦A¥Ñ funcv ¥h©I¥s odeint¡Aodeint ·íµM·|»Ý­n vstart(nvar) ¨Ó§@¶}©l¡A§Ú­Ìª`·N¦b¨ç¦¡ shoot ¡]¤]´N¬O³oùØ newt ¤Uªº funcv¡^¦³¥ý©I¥s¤@­Ó load(x1,v,y) ¡A¦Ó funcv ¥u¦³ V(n2) ¶Ç¶i¨Ó¡A¦]¦¹ ystart(nvar) ¬O­n¦b¥ý§Q¥Î load ¨Ó²Õ¦X¥X¡A½Ò¤å¤¤µ{¦¡»¡©ú¤]¬O³o¼ËÁ¿ªº¡G

x1 »P x2 ¦b¾ã­Ó¹Bºâ¬yµ{¤¤¤@ª½¬O¨D¸Ñ½d³òªº°_©lÂI»P³Ì²×ÂI¨S¦³§ïÅÜ¡Acommon block /caller/ ¦³¥¦©w¸q¡C½Ð¤p¤ß common block ¤¤¦³ªºªF¦è´N¤£¤¹³\¦A¼g¦¨¶Ç¤Þ¼Æªº°Ê§@¤F¡A¤]´N¬O»¡ load »P score ³£¤£¯à¼g¦¨­n¨Ï¥Î /caller/ ¤½¦@°Ï¶ô¡A¦]¬°±q shoot ¥i¬Ý¨£¥¦­Ì¥H¶Ç¤J¤Þ¼Æªº¤è¦¡³B²z x1 »P x2¡A·íµM¨Æ¹ê¤W¤]´N¨S¦³³o­Ó¥²­n common block¡C

¨Ò¦p¡A¨D¸Ñ¤@­Ó¤G¶¥±`·L¤À¤èµ{¦¡¡A¦ÓÃä¬É±ø¥ó¬O y(x1) = 0.0¡By(x2) = 0.0 ®É¡A«h¨ä load °Æµ{¦¡½d¨Ò¦p¤U¡G

subroutine load(x1,v,y)
real x1,v,y(2)
y(1) = 0.0
y(2) = v
end

µù¡G¬°¤°»ò³oùØ­n¶Ç x1 ¶i¥h¡H³o¬O¦]¬°Ãä¬É±ø¥óªº­pºâµøÃD¥Øªº­n¨D¤]¥i¯à·|¥Î¨ì x ­È¡CÁöµM¦b¥»¨Ò¤¤¨S¦³¥Î¨ì¡A¦ý¬°¤F§¹¾ã©Ê¨ä°Æµ{¦¡ªº¤Þ¼Æ¤´­n³]­p¦¨¦³¶Ç¤J¡C

ÁöµM n2 ¦b shoot (funcv) ¤¤¦³©w¸q¡A¦ý«o¨S¦³¥H¤Þ¼Æ¤è¦¡¶Ç¤J score °Æµ{¦¡¡A§A¥i¥H¦Ò¼{ª½±µ§â°}¦C¤j¤p©w¤U¥h©Î¬O¶¡±µ¥t¦æ§ä¤@­Ó·sÅܼƨӨϥΡA¦p¥H¤U½d¨Ò¡]¤@¼ËÃä¬É±ø¥ó¬O y(x1) = 0.0¡By(x2) = 0.0¡^¡G

subroutine score(x2,y,f)
integer n2
parameter (n2=1)
real x2,y(2),f(n2)
f(1) = y(1) - 0.0
write(*,*) 'f value in subroutine score is :', f(1)
end

´£¿ô¤j®a¡A§âÃä¬É±ø¥ó¹F¦¨»P§_ªº«ü¼Ð¨ç¼Æ f ªº­È¦L¥X¨Ó¬O«Ü¦³¥Îªº¡A¥t¥~¡A¦b©I¥s§¹ newt¡]¤]´N¬O§ä§¹¤F®Ú¤§«á¡^¡A§â³Ì²×²Å¦XÃä¬É±ø¥óªº V ­È¦L¥X¤]«Ü¥i¥Hª¾¹D»P²q´ú­È®t¦h¤Ö¡C

¦A¦hµ¹¤@²Õ¨Ò¤l¡A°²³]Ãä¬É±ø¥ó¬° y'(x1) = 1.0¡By'(x2) = 2.5¡A«h load »P score ªº¼gªk¤À§O¬°¡G

subroutine load(x1,v,y)
real x1,v,y(2)
y(1) = v
y(2) = 1.0
end

subroutine score(x2,y,f)
integer n2
parameter (n2=1)
real x2,y(2),f(n2)
f(1) = y(2) - 2.5
write(*,*) 'f value in subroutine score is :', f(1)
end

 

(­×§ï°Æµ{¦¡¤¤ªº³]©w)
ÁöµM½Ò¥»°Æµ{¦¡°ò¥»¤W³£¬O¤£»Ý­n­×§ï´N¥i¥H¨Ï¥Îªº¡A¦ý¥»¸`ªºµ{¦¡²Õ¦X¤F¨â¤j¼Æ­È¤èªk¡A¨ä¤¤¦³«Ü¦h½Õ¾ã©Ê°Ñ¼Æ¦bµ{¦¡½X¤¤³£¬O¥ýµ¹¤F¹w´ú­È¡A¹³¦pªG­n«O¯d¤¤¶¡¹Lµ{ªº¸Ü¡A«e­±ªº´£¨ìªº kmax ´N­n¥´¶}¡]¤]´N¬O»¡¡A¦b shoot.f ¤¤§â­ì¥»ªº kmax = 0 §ï¦¨¬° kmax = KMAXX ¡^¡C¥t¥~¡A±±¨î odeint ºë±K«×ªº EPS ¤]¥iµø§A¦Û¤wªº»Ý­n¦Ó½Õ¾ã¡A¥¦·|ª½±µ¼vÅT¤¤¶¡¹Lµ{ªºÂI¼Æ¡C¦Ü©ó dxsav ¬O½Õ¾ã¦h¤jªº rkqs ¨B´T­È´N­nÀx¦sªº°Ñ¼Æ¡A·Q¦s¦h­ÓÂI´N³]¤p¤@ÂI¡C

(ÅÞ¿èÅÜ¼Æ check)
¦pªGÅÞ¿èÅÜ¼Æ check ªºµ²ªG¬°¯u (T) ¡A¥Nªí­pºâ¦³¨Ç¤£¶¶§Q»Ý­nÀˬd¡F¦pªG¬O°² (F)¡A«h¥¦ªí­pºâ¦³¥¿±`µ²§ô¦Ó¤£»Ý­nÀˬd¡C

(¹ïªì©l­Èªº±Ó·P«×)
¦b¹ê°µªº¹Lµ{¤¤¡A§A·|µo²{¹ï¤£¦Pªº°ÝÃD¹ïªì©l­Èªº±Ó·P«×¬O¤£¤@¼Ëªº¡C¦³¨Ç¤ñ¸û§xÃøªº±¡ªp¡A¤£¾A·íªº²q´úªì©l­È·|¾É­P ludcmp µLªk¨D¸Ñ¡A³q±`¦h¸Õ´X­Óªì©l²q´ú­È¥i¦¸¸Ñ¨M°ÝÃD¡C

(»P PGPLOT ªº·f°t)
­Y­n»P PGPLOT ø¹Ï°Æµ{¦¡®wªº·f°t¡A«h¦b½sĶ®ÉÁÙ­n¦AÃìµ² pgplot »P X11 ¨â­Ó¨ç¦¡®w¡A¨Ò¦p¡G
gfortran -o my_exe.x my_prog.f -lX11 -L /usr/local/pgplot -lpgplot

 

¤@­Ó¹ê»Ú¥i¥Îªº½d¨Ò¥Dµ{¦¡ shoot_main_pgsub.f¡A¦b¦¹¦C¥X¥H¨Ñ°Ñ¦Ò¡G

program shoot_main_pgsub
implicit none
integer i,n2,nvar,kmax,kount,KMAXX,NMAX
parameter (n2=1,NMAX=50,KMAXX=200)
real f(n2),V(n2),x1,x2,dxsav,xp(KMAXX),yp(NMAX,KMAXX)
logical check
common /caller/ x1,x2,nvar
common /path/ kmax,kount,dxsav,xp,yp

nvar = 2
kmax = KMAXX
dxsav = 0.01

write(*,*) 'Please type in beginning and final x : x1 and x2'
read (*,*) x1, x2
write(*,*) 'What is the guess (vector) of the initial condition ?'
read (*,*) (V
(i),i=1,n2)

call newt(V,n2,check)

write(*,*) 'n2 is ;', n2
write(*,*) 'The final value of V is :', V
write(*,*) 'check of newt is', check
write(*,*) 'koumt is', kount

call pg_draw_path()
end

 

subroutine load(x1,V,y)
real x1,V,y(2)
y(1) = 0.0
y(2) = V
end

 

subroutine score(x2,y,f)
integer n2
parameter (n2=1)
real x2,y(2),f(n2)
f(1) = y(1) - 0.0
write(*,*) 'f value in subroutine score is :', f(1)
end

 

subroutine derivs(x,y,dydx)
real x,y(2),dydx(2)
dydx(1) = y(2)
dydx(2) = -y(1)
end

 

subroutine pg_draw_path()
implicit none
integer pgopen,i,j,nvar,nvar1,kmax,kount,KMAXX,NMAX
parameter (nvar1=2,NMAX=50,KMAXX=200)
real x1,x2,dxsav,xp(KMAXX),yp(NMAX,KMAXX)

c For pgplot

real y_min(nvar1),y_max(nvar1),yy(KMAXX),y_extra

common /caller/ x1,x2,nvar
common /path/ kmax,kount,dxsav,xp,yp

c Find Max and min values of yi(x), but here we will use only y1(x).

do i=1,nvar
  y_max(i)= yp(i,1)
  y_min(i)= yp(i,1)
end do

do i=1,nvar
  do j=2,kount
    if (yp(i,j).gt.y_max(i)) y_max(i)=yp(i,j)
    if (yp(i,j).lt.y_min(i)) y_min(i)=yp(i,j)
  end do
end do

c Copy y(1,1:nstep) to yy(1:nstep) for PGLINE plotting

do i=1,kount
  yy(i) = yp(1,i)
end do

c     Intermidiate reult check
c     do i=1,kount
c       write(*,200) i,xp(i),i, yy(i)
c 200   format(1x,'x(',i3,') is ',f5.2,' y(',i3,') is ',f5.2)
c     end do

c Start plotting

y_extra = 0.1*(y_max(1)-y_min(1))
if (pgopen('/xwin').le.0) stop
call pgpap(6.0,0.75)
call pgenv (x1,x2,y_min(1)-y_extra,y_max(1)+y_extra,0,0)
call pgline(kount,xp,yy)
call pgsci(2)
call pgpt(kount,xp,yy,9)
call pgclos

end