#include "rk4.h"

double d1y(double x, double y){
    return pow(x,2) + pow(y,2);
}

double d2y(double x, double y, double z){
    return z;
}

double d2z(double x, double y, double z) {
    return x-3*z-2*y;
}

void rk4_2_variables() {

    double a = 0, b = 1.4, x = 0, y = 0, result = 0, h = 0.1, n = (abs(b-a)/h);
    double d1, d2, d3,d4;

    for(int i = 0; i < n; i++){
        d1 = h*d1y(x, y);
        d2 = h*d1y(x + (h/2), y + (d1/2.0));
        d3 = h*d1y(x +(h/2.0), y + (d2/2.0));
        d4 = h*d1y(x+h, y +d3);
        y += (d1/6.0) + (d2/3.0) + (d3/3.0) + (d4/6.0);
        x += h;
        cout << i << ". x = " << x << "\t y = " << y << endl;
    }

    result = y;

    cout << endl << "---------" << endl;

    x = 0, y = 0;
    double result1, h1 = h/2.0, n1 = abs(b-a)/h1;

    for(int i = 0; i < n1; i++){
        d1 = h1*d1y(x, y);
        d2 = h1*d1y(x + (h1/2), y + (d1/2.0));
        d3 = h1*d1y(x +(h1/2.0), y + (d2/2.0));
        d4 = h1*d1y(x+h1, y +d3);
        y += (d1/6.0) + (d2/3.0) + (d3/3.0) + (d4/6.0);
        x += h1;
        cout << i << ". x = " << x << "\t y = " << y << endl;
    }

    result1 = y;

    cout << endl << "---------" << endl;

    x = 0, y = 0;
    double result2, h2 = h1/2.0, n2 = abs(b-a)/h2;

    for(int i = 0; i < n2; i++){
        d1 = h2*d1y(x, y);
        d2 = h2*d1y(x + (h2/2), y + (d1/2.0));
        d3 = h2*d1y(x +(h2/2.0), y + (d2/2.0));
        d4 = h2*d1y(x+h2, y +d3);
        y += (d1/6.0) + (d2/3.0) + (d3/3.0) + (d4/6.0);
        x += h2;
        cout << i << ". x = " << x << "\t y = " << y << endl;
    }

    result2 = y;

    cout << endl << "---------" << endl;

    double qc = (result1 - result)/(result2 - result1);

    cout << "quociente de convergencia = " << qc << endl;

    double e = result2 - result1;

    cout << "erro = " << e << endl;
}

void rk4_3_variables() {
    double a = 0, b = 0.5, x = 0, y = 0, z = 0, result = 0, h = 0.05, n = (abs(b-a)/h);
    double dy1, dy2, dy3, dy4, dz1, dz2, dz3, dz4;

    for(int i = 0; i < n; i++){
        dy1 = h*d2y(x, y, z);
        dy2 = h*d2y(x + (h/2), y + (dy1/2.0), z + (dy1/2.0));
        dy3 = h*d2y(x +(h/2.0), y + (dy2/2.0), z + (dy2/2.0));
        dy4 = h*d2y(x+h, y +dy3, z + dy3);

        dz1 = h*d2z(x, y, z);
        dz2 = h*d2z(x + (h/2), y + (dz1/2.0), z + (dz1/2.0));
        dz3 = h*d2z(x +(h/2.0), y + (dz2/2.0), z + (dz2/2.0));
        dz4 = h*d2z(x+h, y +dz3, z + dz3);

        y += (dy1/6.0) + (dy2/3.0) + (dy3/3.0) + (dy4/6.0);
        z += (dz1/6.0) + (dz2/3.0) + (dz3/3.0) + (dz4/6.0);
        x += h;


        cout << i << ". x = " << x << "\t y = " << y << "\t z = " << z << endl;
    }
}