Runge Kutta and PGPLOT

131 Views Asked by At

Hi guys when I use this code (makefile not included) I am supposed to be getting two sinusoidal outputs graphed on the same set of axes but for some reason although the output data is accurate the graph is not. In face its not even a smooth curve, it has holes and cusps. Anyone have an idea as to why? Does it have something to do this the array I have written for the time axis? (Note beginner at coding)

#include <iostream> 
#include <cmath>
#include <stdio.h>
#include "cpgplot.h"

//d2xdt2 = -(g/L)*sin(x)
//dwdt = -d2xdt2
//coeff = -(g/L)
//dxdt= w

float dxdt(float w)
        {
            return w;
        }

float dwdt(float x, float coeff)
        {
            return coeff*sin(x);
        }

int main(){

// Standard Variables   
    float wARR[5];  
    float xARR[5]; 
    float tARR[5];    
    float w;       //w = angular velocity
    float x;     //x = theta used in second order equation
    float xo;          //xo = theta initial 
    float wo = 0;       //wo = initial angular velocity 
    float dt = 1;       //h = step value

    float g=9.8;        //g = gravity
    float L;        //L = length of pendulum
    float to = 0;   
    float t;   
    float time;

    float kx1, kw1;
    float kx2, kw2;
    float kx3, kw3;
    float kx4, kw4;

// Input Printing
std::cout << "\n";
std::cout << "Click ENTER key after each value. Please input: \n The initial angle (in decimal radians): "<<"\n";
    std::cin >> xo;
std::cout << "The length of the pendulum (in meters)."<<"\n";
    std::cin >> L;  
    std::cout << "\n";

// Specific Variable Declarations
    float coeff=-(g/L);

//Checkpoint values 
    std::cout << "CHECKPOINT VALUES: ";
    std::cout << "\n";
    std::cout << "Confirming your initial angle is: " << xo << " radian(s)"; 
    std::cout << "\n";
    std::cout << "Confirming the length is: " << L << " meters(s)"; 
    std::cout << "\n";
    std::cout << "Confirming the gravity constant is: " << g << " meters/sec"; 
    std::cout << "\n";
    std::cout << "Coeff (-g/L): " << coeff;
    std::cout << "\n";
    std::cout << "\n";

    std::cout << "Insert the time (in seconds) you want to learn some cool stuff about your pendulum: "<<"\n";
    std::cin >> time;
    std::cout << "\n";

//Array info and Runge Kutta
    xARR[0] = xo;
    wARR[0] =  0;
    tARR[0] = 0;
    time=time+1;

for (int i = 0; i < time ; i++){
    x = xARR[i]; 
    w = wARR[i]; 
    t = tARR[i];

    kx1=dt*dxdt(w);
    kw1=dt*dwdt(x,coeff);
    kx2=dt*dxdt(w+kw1*0.5);
    kw2=dt*dwdt(x+kx1*0.5, coeff);
    kx3=dt*dxdt(w+kw2*0.5);
    kw3=dt*dwdt(x+kx2*0.5, coeff);
    kx4=dt*dxdt(w+kw3);
    kw4=dt*dwdt(x+kx3, coeff);

    std::cout << "\n";  
    std::cout << "RUNGE KUTTA VALUES: at " << i << " second(s)" << "\n";
    std::cout << "kx1:  " << kx1 << "\n";
    std::cout << "kx2:  " << kx2 << "\n";
    std::cout << "kx3:  " << kx3 << "\n";
    std::cout << "kx4:  " << kx4 << "\n";
    std::cout << "kw1:  " << kw1 << "\n";
    std::cout << "kw2:  " << kw2 << "\n";
    std::cout << "kw3:  " << kw3 << "\n";
    std::cout << "kw4:  " << kw4 << "\n";

    xARR[i+1] = xo + (1.0/6.0)*(kx1 + 2*kx2 + 2*kx3 + kx4); 
    wARR[i+1] = wo + (1.0/6.0)*(kw1 + 2*kw2 + 2*kw3 + kw4); 

    std::cout << "FINAL VALUES: at " << i << " second(s)" << "\n";
    std::cout << "The angle is: " << x << " radian(s)"; 
    std::cout << "\n";
    std::cout << "The velocity is: " << w << " radians/s"; 
    std::cout << "\n";
    std::cout << "\n";
    std::cout << "-----------------------------------------------------\n";
    std::cout << "\n";
}
// Graphing with arrays 
//To see the time dependence, let's suppose you have two arrays: one for w and one for x. Then the values x[i] and w[i] entry describe the position and velocity of your pendulum at time t_i.

// Open a plot window
    if (!cpgopen("/XWINDOW")) return 1;

// Set-up plot axes
// M_PI is defined in cmath
    cpgenv(0,time,-1,1,0,1);

// Label axes
    cpglab("time", "angle", "RED = Angle      BLUE = Velocity");

    cpgsci(2); //RED
    cpgline(10,tARR, xARR);
    cpgsci(4); //BLUE
    cpgline(10,tARR, wARR);
    cpgclos();
    }
0

There are 0 best solutions below