## Posted By

Bangonkali on 03/30/12

# 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.

`/*  * File:   main.cpp * Author: Bangonkali ([email protected]) * * Heunâ€™s Method * More indepth discussion on the Book * Numerical Methods 6th Edition - Chapra * Page 722 *  * This example solves the following ordinary differential equation: * y' = ((4 * exp(0.8*x)) - (0.5*y)); 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. *  * Test Run: * at f(4.000000000000000083266726846886740531772375106811523437500000): * estd y: 75.338965077860455168967668271307047689333558082580566406250000 * true y: 75.338962609158578537238426520161738153547048568725585937500000 * true e: 0.000003284801148791693124426737734205834079448393936218053568 * true_max_error: 0.000003284801148791693124426737734205834079448393936218053568 *  * Created on March 29, 2012, 10:03 PM */ #include <cstdlib>#include <iostream>#include <cmath>#include <iomanip> using namespace std; // the formula of y-primelong double f(long double x, long double y) {    return ((4 * exp(0.8*x)) - (0.5*y));} // the correct analytical solutionlong double g(long double x) {    return (4/1.3)*(exp(0.8*x)-exp(-0.5*x))+2*exp(-0.5*x);} // predictorlong double hf0(long double x1, long double y1, long double h) {    return y1 + (f(x1, y1) * h);} // correctorlong double hf0_corrector(long double x1, long double y1, long double x2, long double y2, long double h) {    return y1 + (((f(x1, y1) + f(x2, y2)) / 2)*h);} // get errorlong double get_error(long double a, long double b) {    return abs(((a-b)*100/a));} int main(int argc, char** argv) {    cout << fixed;    cout << setprecision (60);     long double x0 = 0; // initial condition    long double y0 = 2; // initial condition    long double step = 1e-3;  // stepping    long double x_max = 4;  // evaluate integral from x0 to x_max     long double predictor = 0;     long double target_error = 2e-5; // using the error approximation     long double error = 100; // error place holder    long double corrector_new = 0;    unsigned long int i = 0;    long double y_true; // true value of y based on g(x) above    long double e_true = 100; // previous true error    long double e_true_now = 100; // current true error    long double corrector; // corrector place holder     long double true_max_error = 0;     for (x0 = 0; x0<=x_max; x0=x0+step) {        target_error = 2e-5;        error = 100;        e_true = 100;        e_true_now = 101;        corrector_new = 0;        i = 0;         y_true = g(x0 + step);        // cout << "true y: " << y_true << endl;         predictor = hf0(x0, y0, step);        // cout << "predictor = " << predictor << endl;         corrector = hf0_corrector(x0, y0, (x0 + step), predictor, step);        // cout << "corrector = " << corrector << endl;         while (e_true > 1) {            corrector_new = hf0_corrector(x0, y0, (x0 + step), corrector, step);            error = get_error(corrector, corrector_new);            e_true_now = get_error(corrector, y_true);            corrector = corrector_new;            i++;             // cout << "true e: " << e_true << endl;            if (e_true_now == e_true) break;             e_true = e_true_now;        } //        cout << "iterations: " << i << endl;//        cout << "estd y: " << corrector << endl;//        cout << "true e: " << e_true << endl;//        cout << endl;         if (true_max_error < e_true) {            true_max_error = e_true;        }        y0 = corrector;    }     cout << "at f(" << x0 << "): " << endl;    cout << "estd y: " << corrector << endl;    cout << "true y: " << y_true << endl;    cout << "true e: " << e_true << endl;    cout << "true_max_error: " << true_max_error << endl;     return 0;}`