tode.cpp.html - numeric - C++ library with numerical algorithms
(HTM) git clone git://src.adamsgaard.dk/numeric
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tode.cpp.html (17180B)
---
1 <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd">
2 <html>
3 <head>
4 <meta http-equiv="content-type" content="text/html; charset=UTF-8">
5 <title>~/code/numeric/exam/ode.cpp.html</title>
6 <meta name="Generator" content="Vim/7.4">
7 <meta name="plugin-version" content="vim7.4_v1">
8 <meta name="syntax" content="cpp">
9 <meta name="settings" content="number_lines,use_css,pre_wrap,no_foldcolumn,expand_tabs,line_ids,prevent_copy=">
10 <meta name="colorscheme" content="desert">
11 <style type="text/css">
12 <!--
13 pre { white-space: pre-wrap; font-family: monospace; color: #ffffff; background-color: #000000; }
14 body { font-family: monospace; color: #ffffff; background-color: #000000; }
15 * { font-size: 1em; }
16 .Type { color: #008000; }
17 .Statement { color: #804000; }
18 .LineNr { color: #804000; }
19 .PreProc { color: #c000c0; }
20 .Constant { color: #af5f00; }
21 .Special { color: #c000c0; }
22 .Comment { color: #008080; }
23 -->
24 </style>
25
26 <script type='text/javascript'>
27 <!--
28
29 /* function to open any folds containing a jumped-to line before jumping to it */
30 function JumpToLine()
31 {
32 var lineNum;
33 lineNum = window.location.hash;
34 lineNum = lineNum.substr(1); /* strip off '#' */
35
36 if (lineNum.indexOf('L') == -1) {
37 lineNum = 'L'+lineNum;
38 }
39 lineElem = document.getElementById(lineNum);
40 /* Always jump to new location even if the line was hidden inside a fold, or
41 * we corrected the raw number to a line ID.
42 */
43 if (lineElem) {
44 lineElem.scrollIntoView(true);
45 }
46 return true;
47 }
48 if ('onhashchange' in window) {
49 window.onhashchange = JumpToLine;
50 }
51
52 -->
53 </script>
54 </head>
55 <body onload='JumpToLine();'>
56 <pre id='vimCodeElement'>
57 <span id="L1" class="LineNr"> 1 </span><span class="PreProc">#include </span><span class="Constant"><iostream></span>
58 <span id="L2" class="LineNr"> 2 </span><span class="PreProc">#include </span><span class="Constant"><vector></span>
59 <span id="L3" class="LineNr"> 3 </span><span class="PreProc">#include </span><span class="Constant"><cmath></span> <span class="Comment">// for sqrt and pow</span>
60 <span id="L4" class="LineNr"> 4 </span><span class="PreProc">#include </span><span class="Constant"><fstream></span>
61 <span id="L5" class="LineNr"> 5 </span><span class="PreProc">#include </span><span class="Constant">"typedefs.h"</span>
62 <span id="L6" class="LineNr"> 6 </span><span class="PreProc">#include </span><span class="Constant">"ode.h"</span>
63 <span id="L7" class="LineNr"> 7 </span><span class="PreProc">#include </span><span class="Constant">"vector_arithmetic.h"</span>
64 <span id="L8" class="LineNr"> 8 </span>
65 <span id="L9" class="LineNr"> 9 </span><span class="Comment">// Constructor</span>
66 <span id="L10" class="LineNr"> 10 </span>ODE::ODE(std::vector<std::<span class="Type">complex</span><Floattype> >
67 <span id="L11" class="LineNr"> 11 </span> (*f_in)(<span class="Type">const</span> std::<span class="Type">complex</span><Floattype> x,
68 <span id="L12" class="LineNr"> 12 </span> <span class="Type">const</span> std::vector<std::<span class="Type">complex</span><Floattype> > &y),
69 <span id="L13" class="LineNr"> 13 </span> <span class="Type">const</span> std::vector<std::<span class="Type">complex</span><Floattype> > y_start,
70 <span id="L14" class="LineNr"> 14 </span> <span class="Type">const</span> std::<span class="Type">complex</span><Floattype> a_in,
71 <span id="L15" class="LineNr"> 15 </span> <span class="Type">const</span> std::<span class="Type">complex</span><Floattype> b_in,
72 <span id="L16" class="LineNr"> 16 </span> <span class="Type">const</span> Floattype h_start,
73 <span id="L17" class="LineNr"> 17 </span> <span class="Type">const</span> Inttype max_steps,
74 <span id="L18" class="LineNr"> 18 </span> <span class="Type">const</span> Floattype delta_in,
75 <span id="L19" class="LineNr"> 19 </span> <span class="Type">const</span> Floattype epsilon_in,
76 <span id="L20" class="LineNr"> 20 </span> <span class="Type">const</span> Floattype power_in,
77 <span id="L21" class="LineNr"> 21 </span> <span class="Type">const</span> Floattype safety_in)
78 <span id="L22" class="LineNr"> 22 </span> : f(f_in),
79 <span id="L23" class="LineNr"> 23 </span> a(a_in), b(b_in),
80 <span id="L24" class="LineNr"> 24 </span> h((b_in-a_in)*h_start),
81 <span id="L25" class="LineNr"> 25 </span> n_max(max_steps),
82 <span id="L26" class="LineNr"> 26 </span> delta(delta_in), epsilon(epsilon_in),
83 <span id="L27" class="LineNr"> 27 </span> power(power_in), safety(safety_in)
84 <span id="L28" class="LineNr"> 28 </span>{
85 <span id="L29" class="LineNr"> 29 </span> x_list.push_back(a);
86 <span id="L30" class="LineNr"> 30 </span> y_list.push_back(y_start);
87 <span id="L31" class="LineNr"> 31 </span>
88 <span id="L32" class="LineNr"> 32 </span> <span class="Comment">// Perform integration</span>
89 <span id="L33" class="LineNr"> 33 </span> rkdriver();
90 <span id="L34" class="LineNr"> 34 </span>}
91 <span id="L35" class="LineNr"> 35 </span>
92 <span id="L36" class="LineNr"> 36 </span>
93 <span id="L37" class="LineNr"> 37 </span><span class="Comment">// Estimate tolerance</span>
94 <span id="L38" class="LineNr"> 38 </span>Floattype ODE::tau(<span class="Type">const</span> std::vector<std::<span class="Type">complex</span><Floattype> > &y,
95 <span id="L39" class="LineNr"> 39 </span> <span class="Type">const</span> std::<span class="Type">complex</span><Floattype> h_i)
96 <span id="L40" class="LineNr"> 40 </span>{
97 <span id="L41" class="LineNr"> 41 </span> <span class="Statement">return</span> abs((epsilon * cnorm(y) + delta) * sqrt(h_i/(b-a)));
98 <span id="L42" class="LineNr"> 42 </span>}
99 <span id="L43" class="LineNr"> 43 </span>
100 <span id="L44" class="LineNr"> 44 </span><span class="Comment">// Runge-Kutta mid-point stepper</span>
101 <span id="L45" class="LineNr"> 45 </span><span class="Type">void</span> ODE::rkstep12(<span class="Type">const</span> std::<span class="Type">complex</span><Floattype> x0,
102 <span id="L46" class="LineNr"> 46 </span> <span class="Type">const</span> std::vector<std::<span class="Type">complex</span><Floattype> > &y0,
103 <span id="L47" class="LineNr"> 47 </span> std::vector<std::<span class="Type">complex</span><Floattype> > &y1,
104 <span id="L48" class="LineNr"> 48 </span> std::vector<std::<span class="Type">complex</span><Floattype> > &dy)
105 <span id="L49" class="LineNr"> 49 </span>{
106 <span id="L50" class="LineNr"> 50 </span> <span class="Comment">// Denominator 2 used in arithmetic operations</span>
107 <span id="L51" class="LineNr"> 51 </span> <span class="Type">const</span> std::<span class="Type">complex</span><Floattype> den2 (<span class="Constant">2.0f</span>,<span class="Constant">2.0f</span>);
108 <span id="L52" class="LineNr"> 52 </span>
109 <span id="L53" class="LineNr"> 53 </span> <span class="Comment">// Evaluate function at two points</span>
110 <span id="L54" class="LineNr"> 54 </span> (<span class="Type">void</span>)f(x0,y0);
111 <span id="L55" class="LineNr"> 55 </span> <span class="Type">const</span> std::vector<std::<span class="Type">complex</span><Floattype> > k0 = f(x0,y0);
112 <span id="L56" class="LineNr"> 56 </span> <span class="Type">const</span> std::vector<std::<span class="Type">complex</span><Floattype> > k12 = f(x0 + h/den2, y0 + k0*h/den2);
113 <span id="L57" class="LineNr"> 57 </span>
114 <span id="L58" class="LineNr"> 58 </span> <span class="Comment">// Write results to output vectors</span>
115 <span id="L59" class="LineNr"> 59 </span> y1 = y0 + k12*h;
116 <span id="L60" class="LineNr"> 60 </span> dy = (k0 - k12) * h/den2;
117 <span id="L61" class="LineNr"> 61 </span>}
118 <span id="L62" class="LineNr"> 62 </span>
119 <span id="L63" class="LineNr"> 63 </span>
120 <span id="L64" class="LineNr"> 64 </span><span class="Comment">// ODE driver with adaptive step size control</span>
121 <span id="L65" class="LineNr"> 65 </span><span class="Type">void</span> ODE::rkdriver()
122 <span id="L66" class="LineNr"> 66 </span>{
123 <span id="L67" class="LineNr"> 67 </span> std::vector<std::<span class="Type">complex</span><Floattype> > dy(n_max);
124 <span id="L68" class="LineNr"> 68 </span> std::vector<std::<span class="Type">complex</span><Floattype> > y1(n_max);
125 <span id="L69" class="LineNr"> 69 </span>
126 <span id="L70" class="LineNr"> 70 </span> std::<span class="Type">complex</span><Floattype> x;
127 <span id="L71" class="LineNr"> 71 </span> Floattype err, tol;
128 <span id="L72" class="LineNr"> 72 </span> std::vector<std::<span class="Type">complex</span><Floattype> > y;
129 <span id="L73" class="LineNr"> 73 </span>
130 <span id="L74" class="LineNr"> 74 </span> <span class="Statement">while</span> (x_list.back().real() < b.real()
131 <span id="L75" class="LineNr"> 75 </span> || x_list.back().imag() < b.imag())
132 <span id="L76" class="LineNr"> 76 </span> {
133 <span id="L77" class="LineNr"> 77 </span> <span class="Comment">// Get values for this step</span>
134 <span id="L78" class="LineNr"> 78 </span> x = x_list.back();
135 <span id="L79" class="LineNr"> 79 </span> y = y_list.back();
136 <span id="L80" class="LineNr"> 80 </span>
137 <span id="L81" class="LineNr"> 81 </span> <span class="Comment">// Make sure we don't step past the upper boundary</span>
138 <span id="L82" class="LineNr"> 82 </span> <span class="Statement">if</span> ((x + h).real() > b.real()
139 <span id="L83" class="LineNr"> 83 </span> || (x + h).imag() > b.imag())
140 <span id="L84" class="LineNr"> 84 </span> h = b - x;
141 <span id="L85" class="LineNr"> 85 </span>
142 <span id="L86" class="LineNr"> 86 </span> <span class="Comment">// Run Runge-Kutta mid-point stepper</span>
143 <span id="L87" class="LineNr"> 87 </span> rkstep12(x, y, y1, dy);
144 <span id="L88" class="LineNr"> 88 </span>
145 <span id="L89" class="LineNr"> 89 </span> <span class="Comment">// Determine whether the step should be accepted</span>
146 <span id="L90" class="LineNr"> 90 </span> err = cnorm(dy); <span class="Comment">// Error estimate</span>
147 <span id="L91" class="LineNr"> 91 </span> tol = tau(y, h); <span class="Comment">// Tolerance</span>
148 <span id="L92" class="LineNr"> 92 </span> <span class="Statement">if</span> (err < tol) { <span class="Comment">// Step accepted</span>
149 <span id="L93" class="LineNr"> 93 </span> x_list.push_back(x+h);
150 <span id="L94" class="LineNr"> 94 </span> y_list.push_back(y1);
151 <span id="L95" class="LineNr"> 95 </span> }
152 <span id="L96" class="LineNr"> 96 </span>
153 <span id="L97" class="LineNr"> 97 </span> <span class="Comment">// Check that we havn't hit the max. number of steps</span>
154 <span id="L98" class="LineNr"> 98 </span> <span class="Statement">if</span> (x_list.size() == n_max) {
155 <span id="L99" class="LineNr"> 99 </span> std::cerr << <span class="Constant">"Error, the max. number of steps "</span>
156 <span id="L100" class="LineNr">100 </span> << <span class="Constant">"was insufficient</span><span class="Special">\n</span><span class="Constant">"</span>
157 <span id="L101" class="LineNr">101 </span> << <span class="Constant">"Try either increasing the max. number "</span>
158 <span id="L102" class="LineNr">102 </span> << <span class="Constant">"of steps, or decreasing the precision "</span>
159 <span id="L103" class="LineNr">103 </span> << <span class="Constant">"requirements</span><span class="Special">\n</span><span class="Constant">"</span>;
160 <span id="L104" class="LineNr">104 </span> <span class="Statement">return</span>;
161 <span id="L105" class="LineNr">105 </span> }
162 <span id="L106" class="LineNr">106 </span>
163 <span id="L107" class="LineNr">107 </span> <span class="Comment">// Determine new step size</span>
164 <span id="L108" class="LineNr">108 </span> std::<span class="Type">complex</span><Floattype> multiplicator (<span class="Constant">2.0f</span>, <span class="Constant">2.0f</span>);
165 <span id="L109" class="LineNr">109 </span> <span class="Statement">if</span> (err > <span class="Constant">0.0f</span>)
166 <span id="L110" class="LineNr">110 </span> h = h*pow(tol/err, power) * safety;
167 <span id="L111" class="LineNr">111 </span> <span class="Statement">else</span>
168 <span id="L112" class="LineNr">112 </span> h = multiplicator*h;
169 <span id="L113" class="LineNr">113 </span> }
170 <span id="L114" class="LineNr">114 </span>}
171 <span id="L115" class="LineNr">115 </span>
172 <span id="L116" class="LineNr">116 </span>
173 <span id="L117" class="LineNr">117 </span><span class="Comment">// Return the number of steps taken</span>
174 <span id="L118" class="LineNr">118 </span>Inttype ODE::steps()
175 <span id="L119" class="LineNr">119 </span>{
176 <span id="L120" class="LineNr">120 </span> <span class="Statement">return</span> x_list.size();
177 <span id="L121" class="LineNr">121 </span>}
178 <span id="L122" class="LineNr">122 </span>
179 <span id="L123" class="LineNr">123 </span><span class="Type">void</span> ODE::print()
180 <span id="L124" class="LineNr">124 </span>{
181 <span id="L125" class="LineNr">125 </span> <span class="Statement">for</span> (Inttype i=<span class="Constant">0</span>; i<x_list.size(); ++i) {
182 <span id="L126" class="LineNr">126 </span> std::cout << x_list[i] << <span class="Special">'\t'</span>;
183 <span id="L127" class="LineNr">127 </span> <span class="Statement">for</span> (Inttype j=<span class="Constant">0</span>; j<y_list[<span class="Constant">0</span>].size(); ++j)
184 <span id="L128" class="LineNr">128 </span> std::cout << y_list[i][j] << <span class="Special">'\t'</span>;
185 <span id="L129" class="LineNr">129 </span> std::cout << <span class="Special">'\n'</span>;
186 <span id="L130" class="LineNr">130 </span> }
187 <span id="L131" class="LineNr">131 </span>}
188 <span id="L132" class="LineNr">132 </span>
189 <span id="L133" class="LineNr">133 </span><span class="Comment">// Write the x- and y-values to file</span>
190 <span id="L134" class="LineNr">134 </span><span class="Type">void</span> ODE::write(<span class="Type">const</span> <span class="Type">char</span>* filename)
191 <span id="L135" class="LineNr">135 </span>{
192 <span id="L136" class="LineNr">136 </span> std::ofstream outstream;
193 <span id="L137" class="LineNr">137 </span>
194 <span id="L138" class="LineNr">138 </span> <span class="Comment">// Open outfile for write</span>
195 <span id="L139" class="LineNr">139 </span> outstream.open(filename);
196 <span id="L140" class="LineNr">140 </span> <span class="Statement">if</span> (!outstream) {
197 <span id="L141" class="LineNr">141 </span> std::cerr << <span class="Constant">"Error, can't open output file '"</span>
198 <span id="L142" class="LineNr">142 </span> << filename << <span class="Constant">"'.</span><span class="Special">\n</span><span class="Constant">"</span>;
199 <span id="L143" class="LineNr">143 </span> <span class="Statement">return</span>;
200 <span id="L144" class="LineNr">144 </span> }
201 <span id="L145" class="LineNr">145 </span>
202 <span id="L146" class="LineNr">146 </span> <span class="Comment">// Write output values</span>
203 <span id="L147" class="LineNr">147 </span> <span class="Statement">for</span> (Inttype i=<span class="Constant">0</span>; i<x_list.size(); ++i) {
204 <span id="L148" class="LineNr">148 </span> outstream << x_list[i].real() << <span class="Special">'\t'</span>;
205 <span id="L149" class="LineNr">149 </span> outstream << x_list[i].imag() << <span class="Special">'\t'</span>;
206 <span id="L150" class="LineNr">150 </span> <span class="Statement">for</span> (Inttype j=<span class="Constant">0</span>; j<y_list[<span class="Constant">0</span>].size(); ++j) {
207 <span id="L151" class="LineNr">151 </span> outstream << y_list[i][j].real() << <span class="Special">'\t'</span>;
208 <span id="L152" class="LineNr">152 </span> outstream << y_list[i][j].imag() << <span class="Special">'\t'</span>;
209 <span id="L153" class="LineNr">153 </span> }
210 <span id="L154" class="LineNr">154 </span> outstream << <span class="Special">'\n'</span>;
211 <span id="L155" class="LineNr">155 </span> }
212 <span id="L156" class="LineNr">156 </span>
213 <span id="L157" class="LineNr">157 </span> <span class="Comment">// Close file</span>
214 <span id="L158" class="LineNr">158 </span> outstream.close();
215 <span id="L159" class="LineNr">159 </span>
216 <span id="L160" class="LineNr">160 </span> <span class="Statement">if</span> (verbose == <span class="Constant">true</span>)
217 <span id="L161" class="LineNr">161 </span> std::cout << <span class="Constant">"Output written in '"</span> << filename << <span class="Constant">"'.</span><span class="Special">\n</span><span class="Constant">"</span>;
218 <span id="L162" class="LineNr">162 </span>}
219 <span id="L163" class="LineNr">163 </span>
220 </pre>
221 </body>
222 </html>
223 <!-- vim: set foldmethod=manual : -->