Posted By

Bangonkali on 03/30/12


Tagged

method Runge-Kutta Heun


Versions (?)

Runge-Kutta (2nd Order): Heun’s method


 / Published in: C++
 

This example solves the following ordinary differential equation: y' = ((4 * exp(0.8x)) - (0.5y)); graph here: http://j.mp/H3K42D using the Heun’s method. Heun’s method is a second order Runge-Kutta Numerical Method for solving ordinary differential equations. More indepth discussion on the Book Numerical Methods 6th Edition - Chapra Page 722.

  1. /*
  2.  * File: main.cpp
  3.  * Author: Bangonkali ([email protected])
  4.  *
  5.  * Heun’s Method
  6.  * More indepth discussion on the Book
  7.  * Numerical Methods 6th Edition - Chapra
  8.  * Page 722
  9.  *
  10.  * This example solves the following ordinary differential equation:
  11.  * y' = ((4 * exp(0.8*x)) - (0.5*y)); graph here: http://j.mp/H3K42D
  12.  * using the Heun’s method. Heun’s method is a second order Runge-Kutta
  13.  * Numerical Method for solving ordinary differential equations.
  14.  *
  15.  * Test Run:
  16.  * at f(4.000000000000000083266726846886740531772375106811523437500000):
  17.  * estd y: 75.338965077860455168967668271307047689333558082580566406250000
  18.  * true y: 75.338962609158578537238426520161738153547048568725585937500000
  19.  * true e: 0.000003284801148791693124426737734205834079448393936218053568
  20.  * true_max_error: 0.000003284801148791693124426737734205834079448393936218053568
  21.  *
  22.  * Created on March 29, 2012, 10:03 PM
  23.  */
  24.  
  25. #include <cstdlib>
  26. #include <iostream>
  27. #include <cmath>
  28. #include <iomanip>
  29.  
  30. using namespace std;
  31.  
  32. // the formula of y-prime
  33. long double f(long double x, long double y) {
  34. return ((4 * exp(0.8*x)) - (0.5*y));
  35. }
  36.  
  37. // the correct analytical solution
  38. long double g(long double x) {
  39. return (4/1.3)*(exp(0.8*x)-exp(-0.5*x))+2*exp(-0.5*x);
  40. }
  41.  
  42. // predictor
  43. long double hf0(long double x1, long double y1, long double h) {
  44. return y1 + (f(x1, y1) * h);
  45. }
  46.  
  47. // corrector
  48. long double hf0_corrector(long double x1, long double y1, long double x2, long double y2, long double h) {
  49. return y1 + (((f(x1, y1) + f(x2, y2)) / 2)*h);
  50. }
  51.  
  52. // get error
  53. long double get_error(long double a, long double b) {
  54. return abs(((a-b)*100/a));
  55. }
  56.  
  57. int main(int argc, char** argv) {
  58. cout << fixed;
  59. cout << setprecision (60);
  60.  
  61. long double x0 = 0; // initial condition
  62. long double y0 = 2; // initial condition
  63. long double step = 1e-3; // stepping
  64. long double x_max = 4; // evaluate integral from x0 to x_max
  65.  
  66. long double predictor = 0;
  67. long double target_error = 2e-5; // using the error approximation
  68. long double error = 100; // error place holder
  69. long double corrector_new = 0;
  70. unsigned long int i = 0;
  71. long double y_true; // true value of y based on g(x) above
  72. long double e_true = 100; // previous true error
  73. long double e_true_now = 100; // current true error
  74. long double corrector; // corrector place holder
  75.  
  76. long double true_max_error = 0;
  77.  
  78. for (x0 = 0; x0<=x_max; x0=x0+step) {
  79. target_error = 2e-5;
  80. error = 100;
  81. e_true = 100;
  82. e_true_now = 101;
  83. corrector_new = 0;
  84. i = 0;
  85.  
  86. y_true = g(x0 + step);
  87. // cout << "true y: " << y_true << endl;
  88.  
  89. predictor = hf0(x0, y0, step);
  90. // cout << "predictor = " << predictor << endl;
  91.  
  92. corrector = hf0_corrector(x0, y0, (x0 + step), predictor, step);
  93. // cout << "corrector = " << corrector << endl;
  94.  
  95. while (e_true > 1) {
  96. corrector_new = hf0_corrector(x0, y0, (x0 + step), corrector, step);
  97. error = get_error(corrector, corrector_new);
  98. e_true_now = get_error(corrector, y_true);
  99. corrector = corrector_new;
  100. i++;
  101.  
  102. // cout << "true e: " << e_true << endl;
  103. if (e_true_now == e_true) break;
  104.  
  105. e_true = e_true_now;
  106. }
  107.  
  108. // cout << "iterations: " << i << endl;
  109. // cout << "estd y: " << corrector << endl;
  110. // cout << "true e: " << e_true << endl;
  111. // cout << endl;
  112.  
  113. if (true_max_error < e_true) {
  114. true_max_error = e_true;
  115. }
  116. y0 = corrector;
  117. }
  118.  
  119. cout << "at f(" << x0 << "): " << endl;
  120. cout << "estd y: " << corrector << endl;
  121. cout << "true y: " << y_true << endl;
  122. cout << "true e: " << e_true << endl;
  123. cout << "true_max_error: " << true_max_error << endl;
  124.  
  125. return 0;
  126. }

Report this snippet  

You need to login to post a comment.