PageRenderTime 64ms CodeModel.GetById 22ms RepoModel.GetById 0ms app.codeStats 2ms

/doc/html/quadrule_8f90_source.html

http://github.com/ixkael/3DEX
HTML | 6728 lines | 6726 code | 0 blank | 2 comment | 0 complexity | de0964cb4b61e40277ad4dee789592d5 MD5 | raw file
Possible License(s): LGPL-3.0

Large files files are truncated, but you can click here to view the full file

  1. <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
  2. <html xmlns="http://www.w3.org/1999/xhtml">
  3. <head>
  4. <meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
  5. <title>3DEX: /Users/bl/Dropbox/3DEX/src/f90/external/quadrule.f90 Source File</title>
  6. <link href="tabs.css" rel="stylesheet" type="text/css"/>
  7. <link href="search/search.css" rel="stylesheet" type="text/css"/>
  8. <script type="text/javascript" src="search/search.js"></script>
  9. <link href="doxygen.css" rel="stylesheet" type="text/css"/>
  10. </head>
  11. <body onload='searchBox.OnSelectItem(0);'>
  12. <!-- Generated by Doxygen 1.7.4 -->
  13. <script type="text/javascript"><!--
  14. var searchBox = new SearchBox("searchBox", "search",false,'Search');
  15. --></script>
  16. <div id="top">
  17. <div id="titlearea">
  18. <table cellspacing="0" cellpadding="0">
  19. <tbody>
  20. <tr style="height: 56px;">
  21. <td style="padding-left: 0.5em;">
  22. <div id="projectname">3DEX&#160;<span id="projectnumber">1.0</span></div>
  23. <div id="projectbrief">Three-dimensional Fourier-Bessel decomposition</div>
  24. </td>
  25. </tr>
  26. </tbody>
  27. </table>
  28. </div>
  29. <div id="navrow1" class="tabs">
  30. <ul class="tablist">
  31. <li><a href="index.html"><span>Main&#160;Page</span></a></li>
  32. <li><a href="pages.html"><span>Related&#160;Pages</span></a></li>
  33. <li><a href="namespaces.html"><span>Modules</span></a></li>
  34. <li><a href="annotated.html"><span>Data&#160;Types&#160;List</span></a></li>
  35. <li class="current"><a href="files.html"><span>Files</span></a></li>
  36. <li id="searchli">
  37. <div id="MSearchBox" class="MSearchBoxInactive">
  38. <span class="left">
  39. <img id="MSearchSelect" src="search/mag_sel.png"
  40. onmouseover="return searchBox.OnSearchSelectShow()"
  41. onmouseout="return searchBox.OnSearchSelectHide()"
  42. alt=""/>
  43. <input type="text" id="MSearchField" value="Search" accesskey="S"
  44. onfocus="searchBox.OnSearchFieldFocus(true)"
  45. onblur="searchBox.OnSearchFieldFocus(false)"
  46. onkeyup="searchBox.OnSearchFieldChange(event)"/>
  47. </span><span class="right">
  48. <a id="MSearchClose" href="javascript:searchBox.CloseResultsWindow()"><img id="MSearchCloseImg" border="0" src="search/close.png" alt=""/></a>
  49. </span>
  50. </div>
  51. </li>
  52. </ul>
  53. </div>
  54. <div id="navrow2" class="tabs2">
  55. <ul class="tablist">
  56. <li><a href="files.html"><span>File&#160;List</span></a></li>
  57. <li><a href="globals.html"><span>File&#160;Members</span></a></li>
  58. </ul>
  59. </div>
  60. <div class="header">
  61. <div class="headertitle">
  62. <div class="title">/Users/bl/Dropbox/3DEX/src/f90/external/quadrule.f90</div> </div>
  63. </div>
  64. <div class="contents">
  65. <a href="quadrule_8f90.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a><a class="code" href="quadrule_8f90.html#a5db7bba3aa7d37d932b3f545141b52d4">00001</a> <span class="keyword">subroutine </span><a class="code" href="quadrule_8f90.html#a5db7bba3aa7d37d932b3f545141b52d4">bashforth_set</a> ( norder, xtab, weight )
  66. <a name="l00002"></a>00002 <span class="comment">!</span>
  67. <a name="l00003"></a>00003 <span class="comment">!*******************************************************************************</span>
  68. <a name="l00004"></a>00004 <span class="comment">!</span>
  69. <a name="l00005"></a>00005 <span class="comment">!! BASHFORTH_SET sets abscissas and weights for Adams-Bashforth quadrature.</span>
  70. <a name="l00006"></a>00006 <span class="comment">!</span>
  71. <a name="l00007"></a>00007 <span class="comment">!</span>
  72. <a name="l00008"></a>00008 <span class="comment">! Definition:</span>
  73. <a name="l00009"></a>00009 <span class="comment">!</span>
  74. <a name="l00010"></a>00010 <span class="comment">! Adams-Bashforth quadrature formulas are normally used in solving</span>
  75. <a name="l00011"></a>00011 <span class="comment">! ordinary differential equations, and are not really suitable for</span>
  76. <a name="l00012"></a>00012 <span class="comment">! general quadrature computations. However, an Adams-Bashforth formula</span>
  77. <a name="l00013"></a>00013 <span class="comment">! is equivalent to approximating the integral of F(Y(X)) between X(M)</span>
  78. <a name="l00014"></a>00014 <span class="comment">! and X(M+1), using an explicit formula that relies only on known values</span>
  79. <a name="l00015"></a>00015 <span class="comment">! of F(Y(X)) at X(M-N+1) through X(M). For this reason, the formulas</span>
  80. <a name="l00016"></a>00016 <span class="comment">! have been included here.</span>
  81. <a name="l00017"></a>00017 <span class="comment">!</span>
  82. <a name="l00018"></a>00018 <span class="comment">! Suppose the unknown function is denoted by Y(X), with derivative</span>
  83. <a name="l00019"></a>00019 <span class="comment">! F(Y(X)), and that approximate values of the function are known at a</span>
  84. <a name="l00020"></a>00020 <span class="comment">! series of X values, which we write as X(1), X(2), ..., X(M). We write</span>
  85. <a name="l00021"></a>00021 <span class="comment">! the value Y(X(1)) as Y(1) and so on.</span>
  86. <a name="l00022"></a>00022 <span class="comment">!</span>
  87. <a name="l00023"></a>00023 <span class="comment">! Then the solution of the ODE Y&#39;=F(X,Y) at the next point X(M+1) is</span>
  88. <a name="l00024"></a>00024 <span class="comment">! computed by:</span>
  89. <a name="l00025"></a>00025 <span class="comment">!</span>
  90. <a name="l00026"></a>00026 <span class="comment">! Y(M+1) = Y(M) + Integral ( X(M) &lt; X &lt; X(M+1) ) F(Y(X)) dX</span>
  91. <a name="l00027"></a>00027 <span class="comment">! = Y(M) + H * Sum ( 1 &lt;= I &lt;= N ) W(I) * F(Y(M+1-I)) approximately.</span>
  92. <a name="l00028"></a>00028 <span class="comment">!</span>
  93. <a name="l00029"></a>00029 <span class="comment">! In the documentation that follows, we replace F(Y(X)) by F(X).</span>
  94. <a name="l00030"></a>00030 <span class="comment">!</span>
  95. <a name="l00031"></a>00031 <span class="comment">! Integration interval:</span>
  96. <a name="l00032"></a>00032 <span class="comment">!</span>
  97. <a name="l00033"></a>00033 <span class="comment">! [ 0, 1 ].</span>
  98. <a name="l00034"></a>00034 <span class="comment">!</span>
  99. <a name="l00035"></a>00035 <span class="comment">! Weight function:</span>
  100. <a name="l00036"></a>00036 <span class="comment">!</span>
  101. <a name="l00037"></a>00037 <span class="comment">! 1.0D+00</span>
  102. <a name="l00038"></a>00038 <span class="comment">!</span>
  103. <a name="l00039"></a>00039 <span class="comment">! Integral to approximate:</span>
  104. <a name="l00040"></a>00040 <span class="comment">!</span>
  105. <a name="l00041"></a>00041 <span class="comment">! Integral ( 0 &lt;= X &lt;= 1 ) F(X) dX.</span>
  106. <a name="l00042"></a>00042 <span class="comment">!</span>
  107. <a name="l00043"></a>00043 <span class="comment">! Approximate integral:</span>
  108. <a name="l00044"></a>00044 <span class="comment">!</span>
  109. <a name="l00045"></a>00045 <span class="comment">! Sum ( 1 &lt;= I &lt;= NORDER ) WEIGHT(I) * F ( 1 - I ),</span>
  110. <a name="l00046"></a>00046 <span class="comment">!</span>
  111. <a name="l00047"></a>00047 <span class="comment">! Note:</span>
  112. <a name="l00048"></a>00048 <span class="comment">!</span>
  113. <a name="l00049"></a>00049 <span class="comment">! The Adams-Bashforth formulas require equally spaced data.</span>
  114. <a name="l00050"></a>00050 <span class="comment">!</span>
  115. <a name="l00051"></a>00051 <span class="comment">! Here is how the formula is applied in the case with non-unit spacing:</span>
  116. <a name="l00052"></a>00052 <span class="comment">!</span>
  117. <a name="l00053"></a>00053 <span class="comment">! Integral ( A &lt;= X &lt;= A+H ) F(X) dX =</span>
  118. <a name="l00054"></a>00054 <span class="comment">! H * Sum ( 1 &lt;= I &lt;= NORDER ) WEIGHT(I) * F ( A - (I-1)*H ),</span>
  119. <a name="l00055"></a>00055 <span class="comment">! approximately.</span>
  120. <a name="l00056"></a>00056 <span class="comment">!</span>
  121. <a name="l00057"></a>00057 <span class="comment">! The reference lists the second coefficient of the order 8 Adams-Bashforth</span>
  122. <a name="l00058"></a>00058 <span class="comment">! formula as</span>
  123. <a name="l00059"></a>00059 <span class="comment">! weight(2) = -1162169.0D+00 / 120960.0D+00</span>
  124. <a name="l00060"></a>00060 <span class="comment">! but this should be</span>
  125. <a name="l00061"></a>00061 <span class="comment">! weight(2) = -1152169.0D+00 / 120960.0D+00</span>
  126. <a name="l00062"></a>00062 <span class="comment">!</span>
  127. <a name="l00063"></a>00063 <span class="comment">! Reference:</span>
  128. <a name="l00064"></a>00064 <span class="comment">!</span>
  129. <a name="l00065"></a>00065 <span class="comment">! Abramowitz and Stegun,</span>
  130. <a name="l00066"></a>00066 <span class="comment">! Handbook of Mathematical Functions,</span>
  131. <a name="l00067"></a>00067 <span class="comment">! National Bureau of Standards, 1964.</span>
  132. <a name="l00068"></a>00068 <span class="comment">!</span>
  133. <a name="l00069"></a>00069 <span class="comment">! Jean Lapidus and John Seinfeld,</span>
  134. <a name="l00070"></a>00070 <span class="comment">! Numerical Solution of Ordinary Differential Equations,</span>
  135. <a name="l00071"></a>00071 <span class="comment">! Academic Press, 1971.</span>
  136. <a name="l00072"></a>00072 <span class="comment">!</span>
  137. <a name="l00073"></a>00073 <span class="comment">! Daniel Zwillinger, editor,</span>
  138. <a name="l00074"></a>00074 <span class="comment">! Standard Mathematical Tables and Formulae,</span>
  139. <a name="l00075"></a>00075 <span class="comment">! 30th Edition,</span>
  140. <a name="l00076"></a>00076 <span class="comment">! CRC Press, 1996.</span>
  141. <a name="l00077"></a>00077 <span class="comment">!</span>
  142. <a name="l00078"></a>00078 <span class="comment">! Modified:</span>
  143. <a name="l00079"></a>00079 <span class="comment">!</span>
  144. <a name="l00080"></a>00080 <span class="comment">! 15 September 1998</span>
  145. <a name="l00081"></a>00081 <span class="comment">!</span>
  146. <a name="l00082"></a>00082 <span class="comment">! Author:</span>
  147. <a name="l00083"></a>00083 <span class="comment">!</span>
  148. <a name="l00084"></a>00084 <span class="comment">! John Burkardt</span>
  149. <a name="l00085"></a>00085 <span class="comment">!</span>
  150. <a name="l00086"></a>00086 <span class="comment">! Parameters:</span>
  151. <a name="l00087"></a>00087 <span class="comment">!</span>
  152. <a name="l00088"></a>00088 <span class="comment">! Input, integer NORDER, the order of the rule. NORDER should be</span>
  153. <a name="l00089"></a>00089 <span class="comment">! between 1 and 8.</span>
  154. <a name="l00090"></a>00090 <span class="comment">!</span>
  155. <a name="l00091"></a>00091 <span class="comment">! Output, double precision XTAB(NORDER), the abscissas of the rule.</span>
  156. <a name="l00092"></a>00092 <span class="comment">!</span>
  157. <a name="l00093"></a>00093 <span class="comment">! Output, double precision WEIGHT(NORDER), the weights of the rule.</span>
  158. <a name="l00094"></a>00094 <span class="comment">! WEIGHT(1) is the weight at X = 0, WEIGHT(2) the weight at X = -1,</span>
  159. <a name="l00095"></a>00095 <span class="comment">! and so on. The weights are rational, and should sum to 1. Some</span>
  160. <a name="l00096"></a>00096 <span class="comment">! weights may be negative.</span>
  161. <a name="l00097"></a>00097 <span class="comment">!</span>
  162. <a name="l00098"></a>00098 <span class="keyword">implicit none</span>
  163. <a name="l00099"></a>00099 <span class="comment">!</span>
  164. <a name="l00100"></a>00100 <span class="keywordtype">integer</span> norder
  165. <a name="l00101"></a>00101 <span class="comment">!</span>
  166. <a name="l00102"></a>00102 <span class="keywordtype">double precision</span> d
  167. <a name="l00103"></a>00103 <span class="keywordtype">integer</span> i
  168. <a name="l00104"></a>00104 <span class="keywordtype">double precision</span> weight(norder)
  169. <a name="l00105"></a>00105 <span class="keywordtype">double precision</span> xtab(norder)
  170. <a name="l00106"></a>00106 <span class="comment">!</span>
  171. <a name="l00107"></a>00107 <span class="keyword">if</span> ( norder == 1 ) <span class="keyword">then</span>
  172. <a name="l00108"></a>00108
  173. <a name="l00109"></a>00109 weight(1) = 1.0D+00
  174. <a name="l00110"></a>00110
  175. <a name="l00111"></a>00111 <span class="keyword">else</span> <span class="keyword">if</span> ( norder == 2 ) <span class="keyword">then</span>
  176. <a name="l00112"></a>00112
  177. <a name="l00113"></a>00113 d = 2.0D+00
  178. <a name="l00114"></a>00114
  179. <a name="l00115"></a>00115 weight(1) = 3.0D+00 / d
  180. <a name="l00116"></a>00116 weight(2) = - 1.0D+00 / d
  181. <a name="l00117"></a>00117
  182. <a name="l00118"></a>00118 <span class="keyword">else</span> <span class="keyword">if</span> ( norder == 3 ) <span class="keyword">then</span>
  183. <a name="l00119"></a>00119
  184. <a name="l00120"></a>00120 d = 12.0D+00
  185. <a name="l00121"></a>00121
  186. <a name="l00122"></a>00122 weight(1) = 23.0D+00 / d
  187. <a name="l00123"></a>00123 weight(2) = - 16.0D+00 / d
  188. <a name="l00124"></a>00124 weight(3) = 5.0D+00 / d
  189. <a name="l00125"></a>00125
  190. <a name="l00126"></a>00126 <span class="keyword">else</span> <span class="keyword">if</span> ( norder == 4 ) <span class="keyword">then</span>
  191. <a name="l00127"></a>00127
  192. <a name="l00128"></a>00128 d = 24.0D+00
  193. <a name="l00129"></a>00129
  194. <a name="l00130"></a>00130 weight(1) = 55.0D+00 / d
  195. <a name="l00131"></a>00131 weight(2) = - 59.0D+00 / d
  196. <a name="l00132"></a>00132 weight(3) = 37.0D+00 / d
  197. <a name="l00133"></a>00133 weight(4) = - 9.0D+00 / d
  198. <a name="l00134"></a>00134
  199. <a name="l00135"></a>00135 <span class="keyword">else</span> <span class="keyword">if</span> ( norder == 5 ) <span class="keyword">then</span>
  200. <a name="l00136"></a>00136
  201. <a name="l00137"></a>00137 d = 720.0D+00
  202. <a name="l00138"></a>00138
  203. <a name="l00139"></a>00139 weight(1) = 1901.0D+00 / d
  204. <a name="l00140"></a>00140 weight(2) = - 2774.0D+00 / d
  205. <a name="l00141"></a>00141 weight(3) = 2616.0D+00 / d
  206. <a name="l00142"></a>00142 weight(4) = - 1274.0D+00 / d
  207. <a name="l00143"></a>00143 weight(5) = 251.0D+00 / d
  208. <a name="l00144"></a>00144
  209. <a name="l00145"></a>00145 <span class="keyword">else</span> <span class="keyword">if</span> ( norder == 6 ) <span class="keyword">then</span>
  210. <a name="l00146"></a>00146
  211. <a name="l00147"></a>00147 d = 1440.0D+00
  212. <a name="l00148"></a>00148
  213. <a name="l00149"></a>00149 weight(1) = 4277.0D+00 / d
  214. <a name="l00150"></a>00150 weight(2) = - 7923.0D+00 / d
  215. <a name="l00151"></a>00151 weight(3) = 9982.0D+00 / d
  216. <a name="l00152"></a>00152 weight(4) = - 7298.0D+00 / d
  217. <a name="l00153"></a>00153 weight(5) = 2877.0D+00 / d
  218. <a name="l00154"></a>00154 weight(6) = - 475.0D+00 / d
  219. <a name="l00155"></a>00155
  220. <a name="l00156"></a>00156 <span class="keyword">else</span> <span class="keyword">if</span> ( norder == 7 ) <span class="keyword">then</span>
  221. <a name="l00157"></a>00157
  222. <a name="l00158"></a>00158 d = 60480.0D+00
  223. <a name="l00159"></a>00159
  224. <a name="l00160"></a>00160 weight(1) = 198721.0D+00 / d
  225. <a name="l00161"></a>00161 weight(2) = - 447288.0D+00 / d
  226. <a name="l00162"></a>00162 weight(3) = 705549.0D+00 / d
  227. <a name="l00163"></a>00163 weight(4) = - 688256.0D+00 / d
  228. <a name="l00164"></a>00164 weight(5) = 407139.0D+00 / d
  229. <a name="l00165"></a>00165 weight(6) = - 134472.0D+00 / d
  230. <a name="l00166"></a>00166 weight(7) = 19087.0D+00 / d
  231. <a name="l00167"></a>00167
  232. <a name="l00168"></a>00168 <span class="keyword">else</span> <span class="keyword">if</span> ( norder == 8 ) <span class="keyword">then</span>
  233. <a name="l00169"></a>00169
  234. <a name="l00170"></a>00170 d = 120960.0D+00
  235. <a name="l00171"></a>00171
  236. <a name="l00172"></a>00172 weight(1) = 434241.0D+00 / d
  237. <a name="l00173"></a>00173 weight(2) = - 1152169.0D+00 / d
  238. <a name="l00174"></a>00174 weight(3) = 2183877.0D+00 / d
  239. <a name="l00175"></a>00175 weight(4) = - 2664477.0D+00 / d
  240. <a name="l00176"></a>00176 weight(5) = 2102243.0D+00 / d
  241. <a name="l00177"></a>00177 weight(6) = - 1041723.0D+00 / d
  242. <a name="l00178"></a>00178 weight(7) = 295767.0D+00 / d
  243. <a name="l00179"></a>00179 weight(8) = - 36799.0D+00 / d
  244. <a name="l00180"></a>00180
  245. <a name="l00181"></a>00181 <span class="keyword">else</span>
  246. <a name="l00182"></a>00182
  247. <a name="l00183"></a>00183 <span class="keyword">write</span> ( *, <span class="stringliteral">&#39;(a)&#39;</span> ) <span class="stringliteral">&#39; &#39;</span>
  248. <a name="l00184"></a>00184 <span class="keyword">write</span> ( *, <span class="stringliteral">&#39;(a)&#39;</span> ) <span class="stringliteral">&#39;BASHFORTH_SET - Fatal error!&#39;</span>
  249. <a name="l00185"></a>00185 <span class="keyword">write</span> ( *, <span class="stringliteral">&#39;(a,i6)&#39;</span> ) <span class="stringliteral">&#39; Illegal value of NORDER = &#39;</span>, norder
  250. <a name="l00186"></a>00186 <span class="keyword">write</span> ( *, <span class="stringliteral">&#39;(a)&#39;</span> ) <span class="stringliteral">&#39; Legal values are 1 through 8.&#39;</span>
  251. <a name="l00187"></a>00187 stop
  252. <a name="l00188"></a>00188
  253. <a name="l00189"></a>00189 <span class="keyword">end if</span>
  254. <a name="l00190"></a>00190
  255. <a name="l00191"></a>00191 <span class="keyword">do</span> i = 1, norder
  256. <a name="l00192"></a>00192 xtab(i) = dble ( 1 - i )
  257. <a name="l00193"></a>00193 <span class="keyword">end do</span>
  258. <a name="l00194"></a>00194
  259. <a name="l00195"></a>00195 return
  260. <a name="l00196"></a>00196 <span class="keyword">end</span>
  261. <a name="l00197"></a><a class="code" href="quadrule_8f90.html#a33c49c6fa2701ed35a7132dfa2ee7a90">00197</a> <span class="keyword">subroutine </span><a class="code" href="quadrule_8f90.html#a33c49c6fa2701ed35a7132dfa2ee7a90">bdf_set</a> ( norder, alpha, beta, gamma )
  262. <a name="l00198"></a>00198 <span class="comment">!</span>
  263. <a name="l00199"></a>00199 <span class="comment">!*******************************************************************************</span>
  264. <a name="l00200"></a>00200 <span class="comment">!</span>
  265. <a name="l00201"></a>00201 <span class="comment">!! BDF_SET sets weights for backward differentiation ODE weights.</span>
  266. <a name="l00202"></a>00202 <span class="comment">!</span>
  267. <a name="l00203"></a>00203 <span class="comment">!</span>
  268. <a name="l00204"></a>00204 <span class="comment">! Discussion:</span>
  269. <a name="l00205"></a>00205 <span class="comment">!</span>
  270. <a name="l00206"></a>00206 <span class="comment">! GAMMA * Y(N+1) = Sum ( 1 &lt;= I &lt;= NORDER ) ALPHA(I) * Y(N+1-I)</span>
  271. <a name="l00207"></a>00207 <span class="comment">! + dX * BETA * Y&#39;(X(N+1),Y(N+1))</span>
  272. <a name="l00208"></a>00208 <span class="comment">!</span>
  273. <a name="l00209"></a>00209 <span class="comment">! This is equivalent to the backward differentiation corrector formulas.</span>
  274. <a name="l00210"></a>00210 <span class="comment">!</span>
  275. <a name="l00211"></a>00211 <span class="comment">! Modified:</span>
  276. <a name="l00212"></a>00212 <span class="comment">!</span>
  277. <a name="l00213"></a>00213 <span class="comment">! 30 December 1999</span>
  278. <a name="l00214"></a>00214 <span class="comment">!</span>
  279. <a name="l00215"></a>00215 <span class="comment">! Author:</span>
  280. <a name="l00216"></a>00216 <span class="comment">!</span>
  281. <a name="l00217"></a>00217 <span class="comment">! John Burkardt</span>
  282. <a name="l00218"></a>00218 <span class="comment">!</span>
  283. <a name="l00219"></a>00219 <span class="comment">! Parameters:</span>
  284. <a name="l00220"></a>00220 <span class="comment">!</span>
  285. <a name="l00221"></a>00221 <span class="comment">! Input, integer NORDER, the order of the rule, between 1 and 6.</span>
  286. <a name="l00222"></a>00222 <span class="comment">!</span>
  287. <a name="l00223"></a>00223 <span class="comment">! Output, double precision ALPHA(NORDER), BETA, GAMMA, the weights.</span>
  288. <a name="l00224"></a>00224 <span class="comment">!</span>
  289. <a name="l00225"></a>00225 <span class="keyword">implicit none</span>
  290. <a name="l00226"></a>00226 <span class="comment">!</span>
  291. <a name="l00227"></a>00227 <span class="keywordtype">integer</span> norder
  292. <a name="l00228"></a>00228 <span class="comment">!</span>
  293. <a name="l00229"></a>00229 <span class="keywordtype">double precision</span> alpha(norder)
  294. <a name="l00230"></a>00230 <span class="keywordtype">double precision</span> beta
  295. <a name="l00231"></a>00231 <span class="keywordtype">double precision</span> gamma
  296. <a name="l00232"></a>00232 <span class="comment">!</span>
  297. <a name="l00233"></a>00233 <span class="keyword">if</span> ( norder == 1 ) <span class="keyword">then</span>
  298. <a name="l00234"></a>00234 beta = 1.0D+00
  299. <a name="l00235"></a>00235 gamma = 1.0D+00
  300. <a name="l00236"></a>00236 alpha(1) = 1.0D+00
  301. <a name="l00237"></a>00237 <span class="keyword">else</span> <span class="keyword">if</span> ( norder == 2 ) <span class="keyword">then</span>
  302. <a name="l00238"></a>00238 beta = 2.0D+00
  303. <a name="l00239"></a>00239 gamma = 3.0D+00
  304. <a name="l00240"></a>00240 alpha(1) = 4.0D+00
  305. <a name="l00241"></a>00241 alpha(2) = - 1.0D+00
  306. <a name="l00242"></a>00242 <span class="keyword">else</span> <span class="keyword">if</span> ( norder == 3 ) <span class="keyword">then</span>
  307. <a name="l00243"></a>00243 beta = 6.0D+00
  308. <a name="l00244"></a>00244 gamma = 11.0D+00
  309. <a name="l00245"></a>00245 alpha(1) = 18.0D+00
  310. <a name="l00246"></a>00246 alpha(2) = - 9.0D+00
  311. <a name="l00247"></a>00247 alpha(3) = 2.0D+00
  312. <a name="l00248"></a>00248 <span class="keyword">else</span> <span class="keyword">if</span> ( norder == 4 ) <span class="keyword">then</span>
  313. <a name="l00249"></a>00249 beta = 12.0D+00
  314. <a name="l00250"></a>00250 gamma = 25.0D+00
  315. <a name="l00251"></a>00251 alpha(1) = 48.0D+00
  316. <a name="l00252"></a>00252 alpha(2) = - 36.0D+00
  317. <a name="l00253"></a>00253 alpha(3) = 16.0D+00
  318. <a name="l00254"></a>00254 alpha(4) = - 3.0D+00
  319. <a name="l00255"></a>00255 <span class="keyword">else</span> <span class="keyword">if</span> ( norder == 5 ) <span class="keyword">then</span>
  320. <a name="l00256"></a>00256 beta = 60.0D+00
  321. <a name="l00257"></a>00257 gamma = 137.0D+00
  322. <a name="l00258"></a>00258 alpha(1) = 300.0D+00
  323. <a name="l00259"></a>00259 alpha(2) = - 300.0D+00
  324. <a name="l00260"></a>00260 alpha(3) = 200.0D+00
  325. <a name="l00261"></a>00261 alpha(4) = - 75.0D+00
  326. <a name="l00262"></a>00262 alpha(5) = 12.0D+00
  327. <a name="l00263"></a>00263 <span class="keyword">else</span> <span class="keyword">if</span> ( norder == 6 ) <span class="keyword">then</span>
  328. <a name="l00264"></a>00264 beta = 60.0D+00
  329. <a name="l00265"></a>00265 gamma = 147.0D+00
  330. <a name="l00266"></a>00266 alpha(1) = 360.0D+00
  331. <a name="l00267"></a>00267 alpha(2) = - 450.0D+00
  332. <a name="l00268"></a>00268 alpha(3) = 400.0D+00
  333. <a name="l00269"></a>00269 alpha(4) = - 225.0D+00
  334. <a name="l00270"></a>00270 alpha(5) = 72.0D+00
  335. <a name="l00271"></a>00271 alpha(6) = - 10.0D+00
  336. <a name="l00272"></a>00272 <span class="keyword">else</span>
  337. <a name="l00273"></a>00273 <span class="keyword">write</span> ( *, <span class="stringliteral">&#39;(a)&#39;</span> ) <span class="stringliteral">&#39; &#39;</span>
  338. <a name="l00274"></a>00274 <span class="keyword">write</span> ( *, <span class="stringliteral">&#39;(a)&#39;</span> ) <span class="stringliteral">&#39;BDF_SET - Fatal error!&#39;</span>
  339. <a name="l00275"></a>00275 <span class="keyword">write</span> ( *, <span class="stringliteral">&#39;(a,i6)&#39;</span> ) <span class="stringliteral">&#39; Illegal order requested = &#39;</span>, norder
  340. <a name="l00276"></a>00276 stop
  341. <a name="l00277"></a>00277 <span class="keyword">end if</span>
  342. <a name="l00278"></a>00278
  343. <a name="l00279"></a>00279 return
  344. <a name="l00280"></a>00280 <span class="keyword">end</span>
  345. <a name="l00281"></a><a class="code" href="quadrule_8f90.html#a7888ddfbe8865b018471ece2f1ac169b">00281</a> <span class="keyword">subroutine </span><a class="code" href="quadrule_8f90.html#a7888ddfbe8865b018471ece2f1ac169b">bdfc_set</a> ( norder, weight, xtab )
  346. <a name="l00282"></a>00282 <span class="comment">!</span>
  347. <a name="l00283"></a>00283 <span class="comment">!*******************************************************************************</span>
  348. <a name="l00284"></a>00284 <span class="comment">!</span>
  349. <a name="l00285"></a>00285 <span class="comment">!! BDFC_SET sets weights for backward differentiation corrector quadrature.</span>
  350. <a name="l00286"></a>00286 <span class="comment">!</span>
  351. <a name="l00287"></a>00287 <span class="comment">!</span>
  352. <a name="l00288"></a>00288 <span class="comment">! Definition:</span>
  353. <a name="l00289"></a>00289 <span class="comment">!</span>
  354. <a name="l00290"></a>00290 <span class="comment">! A backward differentiation corrector formula is defined for a set</span>
  355. <a name="l00291"></a>00291 <span class="comment">! of evenly spaced abscissas X(I) with X(1) = 1 and X(2) = 0. Assuming</span>
  356. <a name="l00292"></a>00292 <span class="comment">! that the values of the function to be integrated are known at the</span>
  357. <a name="l00293"></a>00293 <span class="comment">! abscissas, the formula is written in terms of the function value at</span>
  358. <a name="l00294"></a>00294 <span class="comment">! X(1), and the backward differences at X(1) that approximate the</span>
  359. <a name="l00295"></a>00295 <span class="comment">! derivatives there.</span>
  360. <a name="l00296"></a>00296 <span class="comment">!</span>
  361. <a name="l00297"></a>00297 <span class="comment">! Integration interval:</span>
  362. <a name="l00298"></a>00298 <span class="comment">!</span>
  363. <a name="l00299"></a>00299 <span class="comment">! [ 0, 1 ]</span>
  364. <a name="l00300"></a>00300 <span class="comment">!</span>
  365. <a name="l00301"></a>00301 <span class="comment">! Weight function:</span>
  366. <a name="l00302"></a>00302 <span class="comment">!</span>
  367. <a name="l00303"></a>00303 <span class="comment">! 1.0D+00</span>
  368. <a name="l00304"></a>00304 <span class="comment">!</span>
  369. <a name="l00305"></a>00305 <span class="comment">! Integral to approximate:</span>
  370. <a name="l00306"></a>00306 <span class="comment">!</span>
  371. <a name="l00307"></a>00307 <span class="comment">! Integral ( 0 &lt;= X &lt;= 1 ) F(X) dX</span>
  372. <a name="l00308"></a>00308 <span class="comment">!</span>
  373. <a name="l00309"></a>00309 <span class="comment">! Approximate integral:</span>
  374. <a name="l00310"></a>00310 <span class="comment">!</span>
  375. <a name="l00311"></a>00311 <span class="comment">! Sum ( 1 &lt;= I &lt;= NORDER ) WEIGHT(I) * BD**(I-1) F ( 1 ).</span>
  376. <a name="l00312"></a>00312 <span class="comment">!</span>
  377. <a name="l00313"></a>00313 <span class="comment">! Here, &quot;BD**(I-1) F ( 1 )&quot; denotes the (I-1)st backward difference</span>
  378. <a name="l00314"></a>00314 <span class="comment">! of F at X = 1, using a spacing of 1. In particular,</span>
  379. <a name="l00315"></a>00315 <span class="comment">!</span>
  380. <a name="l00316"></a>00316 <span class="comment">! BD**0 F(1) = F(1)</span>
  381. <a name="l00317"></a>00317 <span class="comment">! BD**1 F(1) = F(1) - F(0)</span>
  382. <a name="l00318"></a>00318 <span class="comment">! BD**2 F(1) = F(1) - 2 * F(0) + F(-1 )</span>
  383. <a name="l00319"></a>00319 <span class="comment">!</span>
  384. <a name="l00320"></a>00320 <span class="comment">! Note:</span>
  385. <a name="l00321"></a>00321 <span class="comment">!</span>
  386. <a name="l00322"></a>00322 <span class="comment">! The relationship between a backward difference corrector and the</span>
  387. <a name="l00323"></a>00323 <span class="comment">! corresponding Adams-Moulton formula may be illustrated for the</span>
  388. <a name="l00324"></a>00324 <span class="comment">! BDF corrector of order 4:</span>
  389. <a name="l00325"></a>00325 <span class="comment">!</span>
  390. <a name="l00326"></a>00326 <span class="comment">! BD**0 F(1) - 1/2 * BD**1 F(1) - 1/12 * BD**2 F(1) - 1/24 * BDF**3 F(1)</span>
  391. <a name="l00327"></a>00327 <span class="comment">! = F(1)</span>
  392. <a name="l00328"></a>00328 <span class="comment">! - 1/2 * ( F(1) - F(0) )</span>
  393. <a name="l00329"></a>00329 <span class="comment">! - 1/12 * ( F(1) - 2 * F(0) + F(-1) )</span>
  394. <a name="l00330"></a>00330 <span class="comment">! - 1/24 * ( F(1) - 3 * F(0) + 3 * F(-1) - F(-2) )</span>
  395. <a name="l00331"></a>00331 <span class="comment">! = 9/24 * F(1) + 19/24 * F(0) - 5/24 * F(-1) + 1/24 * F(-2)</span>
  396. <a name="l00332"></a>00332 <span class="comment">!</span>
  397. <a name="l00333"></a>00333 <span class="comment">! which is the Adams-Moulton formula of order 4.</span>
  398. <a name="l00334"></a>00334 <span class="comment">! </span>
  399. <a name="l00335"></a>00335 <span class="comment">! Reference:</span>
  400. <a name="l00336"></a>00336 <span class="comment">!</span>
  401. <a name="l00337"></a>00337 <span class="comment">! Simeon Fatunla,</span>
  402. <a name="l00338"></a>00338 <span class="comment">! Numerical Methods for Initial Value Problems in Ordinary Differential</span>
  403. <a name="l00339"></a>00339 <span class="comment">! Equations,</span>
  404. <a name="l00340"></a>00340 <span class="comment">! Academic Press, 1988.</span>
  405. <a name="l00341"></a>00341 <span class="comment">!</span>
  406. <a name="l00342"></a>00342 <span class="comment">! Modified:</span>
  407. <a name="l00343"></a>00343 <span class="comment">!</span>
  408. <a name="l00344"></a>00344 <span class="comment">! 28 February 2000</span>
  409. <a name="l00345"></a>00345 <span class="comment">!</span>
  410. <a name="l00346"></a>00346 <span class="comment">! Author:</span>
  411. <a name="l00347"></a>00347 <span class="comment">!</span>
  412. <a name="l00348"></a>00348 <span class="comment">! John Burkardt</span>
  413. <a name="l00349"></a>00349 <span class="comment">!</span>
  414. <a name="l00350"></a>00350 <span class="comment">! Parameters:</span>
  415. <a name="l00351"></a>00351 <span class="comment">!</span>
  416. <a name="l00352"></a>00352 <span class="comment">! Input, integer NORDER, the order of the rule, which can be</span>
  417. <a name="l00353"></a>00353 <span class="comment">! any value from 1 to 19.</span>
  418. <a name="l00354"></a>00354 <span class="comment">!</span>
  419. <a name="l00355"></a>00355 <span class="comment">! Output, double precision WEIGHT(NORDER), the weights of the rule.</span>
  420. <a name="l00356"></a>00356 <span class="comment">!</span>
  421. <a name="l00357"></a>00357 <span class="comment">! Output, double precision XTAB(NORDER), the abscissas of the rule.</span>
  422. <a name="l00358"></a>00358 <span class="comment">!</span>
  423. <a name="l00359"></a>00359 <span class="keyword">implicit none</span>
  424. <a name="l00360"></a>00360 <span class="comment">!</span>
  425. <a name="l00361"></a>00361 <span class="keywordtype">integer</span>, <span class="keywordtype">parameter</span> :: maxord = 19
  426. <a name="l00362"></a>00362 <span class="comment">!</span>
  427. <a name="l00363"></a>00363 <span class="keywordtype">integer</span> norder
  428. <a name="l00364"></a>00364 <span class="comment">!</span>
  429. <a name="l00365"></a>00365 <span class="keywordtype">integer</span> i
  430. <a name="l00366"></a>00366 <span class="keywordtype">double precision</span> w(maxord)
  431. <a name="l00367"></a>00367 <span class="keywordtype">double precision</span> weight(norder)
  432. <a name="l00368"></a>00368 <span class="keywordtype">double precision</span> xtab(norder)
  433. <a name="l00369"></a>00369 <span class="comment">!</span>
  434. <a name="l00370"></a>00370 w(1) = 1.0D+00
  435. <a name="l00371"></a>00371 w(2) = - 1.0D+00 / 2.0D+00
  436. <a name="l00372"></a>00372 w(3) = - 1.0D+00 / 12.0D+00
  437. <a name="l00373"></a>00373 w(4) = - 1.0D+00 / 24.0D+00
  438. <a name="l00374"></a>00374 w(5) = - 19.0D+00 / 720.0D+00
  439. <a name="l00375"></a>00375 w(6) = - 3.0D+00 / 160.0D+00
  440. <a name="l00376"></a>00376 w(7) = - 863.0D+00 / 60480.0D+00
  441. <a name="l00377"></a>00377 w(8) = - 275.0D+00 / 24792.0D+00
  442. <a name="l00378"></a>00378 w(9) = - 33953.0D+00 / 3628800.0D+00
  443. <a name="l00379"></a>00379 w(10) = - 8183.0D+00 / 1036800.0D+00
  444. <a name="l00380"></a>00380 w(11) = - 3250433.0D+00 / 479001600.0D+00
  445. <a name="l00381"></a>00381 w(12) = - 4671.0D+00 / 788480.0D+00
  446. <a name="l00382"></a>00382 w(13) = - 13695779093.0D+00 / 2615348736000.0D+00
  447. <a name="l00383"></a>00383 w(14) = - 2224234463.0D+00 / 475517952000.0D+00
  448. <a name="l00384"></a>00384 w(15) = - 132282840127.0D+00 / 31384184832000.0D+00
  449. <a name="l00385"></a>00385 w(16) = - 2639651053.0D+00 / 689762304000.0D+00
  450. <a name="l00386"></a>00386 w(17) = 111956703448001.0D+00 / 3201186852864.0D+00
  451. <a name="l00387"></a>00387 w(18) = 50188465.0D+00 / 15613165568.0D+00
  452. <a name="l00388"></a>00388 w(19) = 2334028946344463.0D+00 / 786014494949376.0D+00
  453. <a name="l00389"></a>00389
  454. <a name="l00390"></a>00390 <span class="keyword">do</span> i = 1, min ( norder, maxord )
  455. <a name="l00391"></a>00391 weight(i) = w(i)
  456. <a name="l00392"></a>00392 <span class="keyword">end do</span>
  457. <a name="l00393"></a>00393
  458. <a name="l00394"></a>00394 <span class="keyword">do</span> i = 1, norder
  459. <a name="l00395"></a>00395 xtab(i) = dble ( 2 - i )
  460. <a name="l00396"></a>00396 <span class="keyword">end do</span>
  461. <a name="l00397"></a>00397
  462. <a name="l00398"></a>00398 return
  463. <a name="l00399"></a>00399 <span class="keyword">end</span>
  464. <a name="l00400"></a><a class="code" href="quadrule_8f90.html#ae1f2722140b0af5e1a2bf8e5877222aa">00400</a> <span class="keyword">subroutine </span><a class="code" href="quadrule_8f90.html#ae1f2722140b0af5e1a2bf8e5877222aa">bdfp_set</a> ( norder, weight, xtab )
  465. <a name="l00401"></a>00401 <span class="comment">!</span>
  466. <a name="l00402"></a>00402 <span class="comment">!*******************************************************************************</span>
  467. <a name="l00403"></a>00403 <span class="comment">!</span>
  468. <a name="l00404"></a>00404 <span class="comment">!! BDFP_SET sets weights for backward differentiation predictor quadrature.</span>
  469. <a name="l00405"></a>00405 <span class="comment">!</span>
  470. <a name="l00406"></a>00406 <span class="comment">!</span>
  471. <a name="l00407"></a>00407 <span class="comment">! Definition:</span>
  472. <a name="l00408"></a>00408 <span class="comment">!</span>
  473. <a name="l00409"></a>00409 <span class="comment">! A backward differentiation predictor formula is defined for a set</span>
  474. <a name="l00410"></a>00410 <span class="comment">! of evenly spaced abscissas X(I) with X(1) = 1 and X(2) = 0. Assuming</span>
  475. <a name="l00411"></a>00411 <span class="comment">! that the values of the function to be integrated are known at the</span>
  476. <a name="l00412"></a>00412 <span class="comment">! abscissas, the formula is written in terms of the function value at</span>
  477. <a name="l00413"></a>00413 <span class="comment">! X(2), and the backward differences at X(2) that approximate the</span>
  478. <a name="l00414"></a>00414 <span class="comment">! derivatives there. A backward differentiation predictor formula</span>
  479. <a name="l00415"></a>00415 <span class="comment">! is equivalent to an Adams-Bashforth formula of the same order.</span>
  480. <a name="l00416"></a>00416 <span class="comment">!</span>
  481. <a name="l00417"></a>00417 <span class="comment">! Integration interval:</span>
  482. <a name="l00418"></a>00418 <span class="comment">!</span>
  483. <a name="l00419"></a>00419 <span class="comment">! [ 0, 1 ]</span>
  484. <a name="l00420"></a>00420 <span class="comment">!</span>
  485. <a name="l00421"></a>00421 <span class="comment">! Weight function:</span>
  486. <a name="l00422"></a>00422 <span class="comment">!</span>
  487. <a name="l00423"></a>00423 <span class="comment">! 1.0D+00</span>
  488. <a name="l00424"></a>00424 <span class="comment">!</span>
  489. <a name="l00425"></a>00425 <span class="comment">! Integral to approximate:</span>
  490. <a name="l00426"></a>00426 <span class="comment">!</span>
  491. <a name="l00427"></a>00427 <span class="comment">! Integral ( 0 &lt;= X &lt;= 1 ) F(X) dX</span>
  492. <a name="l00428"></a>00428 <span class="comment">!</span>
  493. <a name="l00429"></a>00429 <span class="comment">! Approximate integral:</span>
  494. <a name="l00430"></a>00430 <span class="comment">!</span>
  495. <a name="l00431"></a>00431 <span class="comment">! Sum ( 1 &lt;= I &lt;= NORDER ) WEIGHT(I) * BD**(I-1) F ( 0 ),</span>
  496. <a name="l00432"></a>00432 <span class="comment">!</span>
  497. <a name="l00433"></a>00433 <span class="comment">! Here, &quot;BD**(I-1) F ( 0 )&quot; denotes the (I-1)st backward difference</span>
  498. <a name="l00434"></a>00434 <span class="comment">! of F at X = 0, using a spacing of 1. In particular,</span>
  499. <a name="l00435"></a>00435 <span class="comment">!</span>
  500. <a name="l00436"></a>00436 <span class="comment">! BD**0 F(0) = F(0)</span>
  501. <a name="l00437"></a>00437 <span class="comment">! BD**1 F(0) = F(0) - F(-1)</span>
  502. <a name="l00438"></a>00438 <span class="comment">! BD**2 F(0) = F(0) - 2 * F(-1) + F(-2 )</span>
  503. <a name="l00439"></a>00439 <span class="comment">!</span>
  504. <a name="l00440"></a>00440 <span class="comment">! Note:</span>
  505. <a name="l00441"></a>00441 <span class="comment">!</span>
  506. <a name="l00442"></a>00442 <span class="comment">! The relationship between a backward difference predictor and the</span>
  507. <a name="l00443"></a>00443 <span class="comment">! corresponding Adams-Bashforth formula may be illustrated for the</span>
  508. <a name="l00444"></a>00444 <span class="comment">! BDF predictor of order 3:</span>
  509. <a name="l00445"></a>00445 <span class="comment">!</span>
  510. <a name="l00446"></a>00446 <span class="comment">! BD**0 F(0) + 0.5 * BD**1 F(0) + 5/12 * BD**2 F(0)</span>
  511. <a name="l00447"></a>00447 <span class="comment">! = F(0)</span>
  512. <a name="l00448"></a>00448 <span class="comment">! + 1/2 * ( F(0) - F(1) )</span>
  513. <a name="l00449"></a>00449 <span class="comment">! + 5/12 * ( F(0) - 2 * F(-1) + F(-2) )</span>
  514. <a name="l00450"></a>00450 <span class="comment">! = 23/12 * F(0) - 16/12 * F(-1) + 5/12 F(-2)</span>
  515. <a name="l00451"></a>00451 <span class="comment">!</span>
  516. <a name="l00452"></a>00452 <span class="comment">! which is the Adams-Bashforth formula of order 3.</span>
  517. <a name="l00453"></a>00453 <span class="comment">! </span>
  518. <a name="l00454"></a>00454 <span class="comment">! Reference:</span>
  519. <a name="l00455"></a>00455 <span class="comment">!</span>
  520. <a name="l00456"></a>00456 <span class="comment">! Simeon Fatunla,</span>
  521. <a name="l00457"></a>00457 <span class="comment">! Numerical Methods for Initial Value Problems in Ordinary Differential</span>
  522. <a name="l00458"></a>00458 <span class="comment">! Equations,</span>
  523. <a name="l00459"></a>00459 <span class="comment">! Academic Press, 1988.</span>
  524. <a name="l00460"></a>00460 <span class="comment">!</span>
  525. <a name="l00461"></a>00461 <span class="comment">! Modified:</span>
  526. <a name="l00462"></a>00462 <span class="comment">!</span>
  527. <a name="l00463"></a>00463 <span class="comment">! 29 February 2000</span>
  528. <a name="l00464"></a>00464 <span class="comment">!</span>
  529. <a name="l00465"></a>00465 <span class="comment">! Author:</span>
  530. <a name="l00466"></a>00466 <span class="comment">!</span>
  531. <a name="l00467"></a>00467 <span class="comment">! John Burkardt</span>
  532. <a name="l00468"></a>00468 <span class="comment">!</span>
  533. <a name="l00469"></a>00469 <span class="comment">! Parameters:</span>
  534. <a name="l00470"></a>00470 <span class="comment">!</span>
  535. <a name="l00471"></a>00471 <span class="comment">! Input, integer NORDER, the order of the rule, which can be</span>
  536. <a name="l00472"></a>00472 <span class="comment">! any value from 1 to 19.</span>
  537. <a name="l00473"></a>00473 <span class="comment">!</span>
  538. <a name="l00474"></a>00474 <span class="comment">! Output, double precision WEIGHT(NORDER), the weight of the rule.</span>
  539. <a name="l00475"></a>00475 <span class="comment">!</span>
  540. <a name="l00476"></a>00476 <span class="comment">! Output, double precision XTAB(NORDER), the abscissas of the rule.</span>
  541. <a name="l00477"></a>00477 <span class="comment">!</span>
  542. <a name="l00478"></a>00478 <span class="keyword">implicit none</span>
  543. <a name="l00479"></a>00479 <span class="comment">!</span>
  544. <a name="l00480"></a>00480 <span class="keywordtype">integer</span>, <span class="keywordtype">parameter</span> :: maxord = 19
  545. <a name="l00481"></a>00481 <span class="comment">!</span>
  546. <a name="l00482"></a>00482 <span class="keywordtype">integer</span> norder
  547. <a name="l00483"></a>00483 <span class="comment">!</span>
  548. <a name="l00484"></a>00484 <span class="keywordtype">integer</span> i
  549. <a name="l00485"></a>00485 <span class="keywordtype">double precision</span> w(maxord)
  550. <a name="l00486"></a>00486 <span class="keywordtype">double precision</span> weight(norder)
  551. <a name="l00487"></a>00487 <span class="keywordtype">double precision</span> xtab(norder)
  552. <a name="l00488"></a>00488 <span class="comment">!</span>
  553. <a name="l00489"></a>00489 w(1) = 1.0D+00
  554. <a name="l00490"></a>00490 w(2) = 1.0D+00 / 2.0D+00
  555. <a name="l00491"></a>00491 w(3) = 5.0D+00 / 12.0D+00
  556. <a name="l00492"></a>00492 w(4) = 3.0D+00 / 8.0D+00
  557. <a name="l00493"></a>00493 w(5) = 251.0D+00 / 720.0D+00
  558. <a name="l00494"></a>00494 w(6) = 95.0D+00 / 288.0D+00
  559. <a name="l00495"></a>00495 w(7) = 19087.0D+00 / 60480.0D+00
  560. <a name="l00496"></a>00496 w(8) = 5257.0D+00 / 17280.0D+00
  561. <a name="l00497"></a>00497 w(9) = 1070017.0D+00 / 3628800.0D+00
  562. <a name="l00498"></a>00498 w(10) = 25713.0D+00 / 89600.0D+00
  563. <a name="l00499"></a>00499 w(11) = 26842253.0D+00 / 95800320.0D+00
  564. <a name="l00500"></a>00500 w(12) = 4777223.0D+00 / 17418240.0D+00
  565. <a name="l00501"></a>00501 w(13) = 703604254357.0D+00 / 2615348736000.0D+00
  566. <a name="l00502"></a>00502 w(14) = 106364763817.0D+00 / 402361344000.0D+00
  567. <a name="l00503"></a>00503 w(15) = 1166309819657.0D+00 / 4483454976000.0D+00
  568. <a name="l00504"></a>00504 w(16) = 25221445.0D+00 / 98402304.0D+00
  569. <a name="l00505"></a>00505 w(17) = 8092989203533249.0D+00 / 3201186852864.0D+00
  570. <a name="l00506"></a>00506 w(18) = 85455477715379.0D+00 / 34237292544.0D+00
  571. <a name="l00507"></a>00507 w(19) = 12600467236042756559.0D+00 / 5109094217170944.0D+00
  572. <a name="l00508"></a>00508
  573. <a name="l00509"></a>00509 <span class="keyword">do</span> i = 1, min ( norder, maxord )
  574. <a name="l00510"></a>00510 weight(i) = w(i)
  575. <a name="l00511"></a>00511 <span class="keyword">end do</span>
  576. <a name="l00512"></a>00512
  577. <a name="l00513"></a>00513 <span class="keyword">do</span> i = 1, norder
  578. <a name="l00514"></a>00514 xtab(i) = dble ( 1 - i )
  579. <a name="l00515"></a>00515 <span class="keyword">end do</span>
  580. <a name="l00516"></a>00516
  581. <a name="l00517"></a>00517 return
  582. <a name="l00518"></a>00518 <span class="keyword">end</span>
  583. <a name="l00519"></a><a class="code" href="quadrule_8f90.html#aed98218103418a57a42bce561d3f221a">00519</a> <span class="keyword">subroutine </span><a class="code" href="quadrule_8f90.html#aed98218103418a57a42bce561d3f221a">bdf_sum</a> ( func, norder, weight, xtab, diftab, result )
  584. <a name="l00520"></a>00520 <span class="comment">!</span>
  585. <a name="l00521"></a>00521 <span class="comment">!*******************************************************************************</span>
  586. <a name="l00522"></a>00522 <span class="comment">!</span>
  587. <a name="l00523"></a>00523 <span class="comment">!! BDF_SUM carries out an explicit backward difference quadrature rule for [0,1].</span>
  588. <a name="l00524"></a>00524 <span class="comment">!</span>
  589. <a name="l00525"></a>00525 <span class="comment">!</span>
  590. <a name="l00526"></a>00526 <span class="comment">! Integral to approximate:</span>
  591. <a name="l00527"></a>00527 <span class="comment">!</span>
  592. <a name="l00528"></a>00528 <span class="comment">! Integral ( 0 &lt;= X &lt;= 1 ) F(X) dX</span>
  593. <a name="l00529"></a>00529 <span class="comment">!</span>
  594. <a name="l00530"></a>00530 <span class="comment">! Formula:</span>
  595. <a name="l00531"></a>00531 <span class="comment">!</span>
  596. <a name="l00532"></a>00532 <span class="comment">! RESULT = Sum ( 1 &lt;= I &lt;= NORDER ) WEIGHT(I) * BDF**(I-1) FUNC ( 0 )</span>
  597. <a name="l00533"></a>00533 <span class="comment">!</span>
  598. <a name="l00534"></a>00534 <span class="comment">! Note:</span>
  599. <a name="l00535"></a>00535 <span class="comment">!</span>
  600. <a name="l00536"></a>00536 <span class="comment">! The integral from 0 to 1 is approximated using data at X = 0,</span>
  601. <a name="l00537"></a>00537 <span class="comment">! -1, -2, ..., -NORDER+1. This is a form of extrapolation, and</span>
  602. <a name="l00538"></a>00538 <span class="comment">! the approximation can become poor as NORDER increases.</span>
  603. <a name="l00539"></a>00539 <span class="comment">!</span>
  604. <a name="l00540"></a>00540 <span class="comment">! Modified:</span>
  605. <a name="l00541"></a>00541 <span class="comment">!</span>
  606. <a name="l00542"></a>00542 <span class="comment">! 26 October 2000</span>
  607. <a name="l00543"></a>00543 <span class="comment">!</span>
  608. <a name="l00544"></a>00544 <span class="comment">! Author:</span>
  609. <a name="l00545"></a>00545 <span class="comment">!</span>
  610. <a name="l00546"></a>00546 <span class="comment">! John Burkardt</span>
  611. <a name="l00547"></a>00547 <span class="comment">!</span>
  612. <a name="l00548"></a>00548 <span class="comment">! Parameters:</span>
  613. <a name="l00549"></a>00549 <span class="comment">!</span>
  614. <a name="l00550"></a>00550 <span class="comment">! Input, external FUNC, the name of the FORTRAN function which evaluates</span>
  615. <a name="l00551"></a>00551 <span class="comment">! the integrand. The function must have the form</span>
  616. <a name="l00552"></a>00552 <span class="comment">! double precision func ( x ).</span>
  617. <a name="l00553"></a>00553 <span class="comment">!</span>
  618. <a name="l00554"></a>00554 <span class="comment">! Input, integer NORDER, the order of the rule.</span>
  619. <a name="l00555"></a>00555 <span class="comment">!</span>
  620. <a name="l00556"></a>00556 <span class="comment">! Input, double precision WEIGHT(NORDER), the weights of the rule.</span>
  621. <a name="l00557"></a>00557 <span class="comment">!</span>
  622. <a name="l00558"></a>00558 <span class="comment">! Input, double precision XTAB(NORDER), the abscissas of the rule.</span>
  623. <a name="l00559"></a>00559 <span class="comment">!</span>
  624. <a name="l00560"></a>00560 <span class="comment">! Workspace, double precision DIFTAB(NORDER).</span>
  625. <a name="l00561"></a>00561 <span class="comment">!</span>
  626. <a name="l00562"></a>00562 <span class="comment">! Output, double precision RESULT, the approximate value of the integral.</span>
  627. <a name="l00563"></a>00563 <span class="comment">!</span>
  628. <a name="l00564"></a>00564 <span class="keyword">implicit none</span>
  629. <a name="l00565"></a>00565 <span class="comment">!</span>
  630. <a name="l00566"></a>00566 <span class="keywordtype">integer</span> norder
  631. <a name="l00567"></a>00567 <span class="comment">!</span>
  632. <a name="l00568"></a>00568 <span class="keywordtype">double precision</span> diftab(norder)
  633. <a name="l00569"></a>00569 <span class="keywordtype">double precision</span>, <span class="keywordtype">external</span> :: func
  634. <a name="l00570"></a>00570 <span class="keywordtype">integer</span> i
  635. <a name="l00571"></a>00571 <span class="keywordtype">integer</span> j
  636. <a name="l00572"></a>00572 <span class="keywordtype">double precision</span> result
  637. <a name="l00573"></a>00573 <span class="keywordtype">double precision</span> weight(norder)
  638. <a name="l00574"></a>00574 <span class="keywordtype">double precision</span> xtab(norder)
  639. <a name="l00575"></a>00575 <span class="comment">!</span>
  640. <a name="l00576"></a>00576 <span class="keyword">do</span> i = 1, norder
  641. <a name="l00577"></a>00577 diftab(i) = func ( xtab(i) )
  642. <a name="l00578"></a>00578 <span class="keyword">end do</span>
  643. <a name="l00579"></a>00579
  644. <a name="l00580"></a>00580 <span class="keyword">do</span> i = 2, norder
  645. <a name="l00581"></a>00581 <span class="keyword">do</span> j = i, norder
  646. <a name="l00582"></a>00582 diftab(norder+i-j) = ( diftab(norder+i-j-1) - diftab(norder+i-j) )
  647. <a name="l00583"></a>00583 <span class="keyword">end do</span>
  648. <a name="l00584"></a>00584 <span class="keyword">end do</span>
  649. <a name="l00585"></a>00585
  650. <a name="l00586"></a>00586 result = dot_product ( weight(1:norder), diftab(1:norder) )
  651. <a name="l00587"></a>00587
  652. <a name="l00588"></a>00588 return
  653. <a name="l00589"></a>00589 <span class="keyword">end</span>
  654. <a name="l00590"></a><a class="code" href="quadrule_8f90.html#ade9f4e674221ccedd6251da05db6bbba">00590</a> <span class="keyword">subroutine </span><a class="code" href="quadrule_8f90.html#ade9f4e674221ccedd6251da05db6bbba">cheb_set</a> ( norder, xtab, weight )
  655. <a name="l00591"></a>00591 <span class="comment">!</span>
  656. <a name="l00592"></a>00592 <span class="comment">!*******************************************************************************</span>
  657. <a name="l00593"></a>00593 <span class="comment">!</span>
  658. <a name="l00594"></a>00594 <span class="comment">!! CHEB_SET sets abscissas and weights for Chebyshev quadrature.</s

Large files files are truncated, but you can click here to view the full file