0
$\begingroup$

I am trying to solve equations of motion for an harmonic oscillator using 4th order runge kutta method, but as a result I get almost constant velocity and position; I feel that the problem is that I did not fully understood the method.

This is th code I am using (in c++)

#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <fstream>
#include <iomanip>

using namespace std;

ofstream output;

void rk(double,double*,double*, double*,double*, int);
void print_file(double*, double*,double*,int n);

int main(int argc, char* argv[]){
    double *t,*x,*v,*ic;

    double h,in=0,fin=2*M_PI;
    int n=10000;

    char* file=argv[1];
    output.open(file);

    h=(fin-in)/n;

    t=new double[n];
    x=new double[n];
    v=new double[n];
    ic=new double[2];

    ic[0]=1;//in. cond. pos.
    ic[1]=0;//in. cond. vel.

    rk(fun,h,ic,t,v,x,n);
    for(int i=0;i<n;i+=500){cout << x[i] << "  " << v[i] << endl;}

    print_file(t,x,v,n);

    delete[] t;
    delete[] x;
    delete[] v;
    delete[] ic;
}

void rk(double step,double *ic,double* t, double* v, double* x, int n){
    double k1,k2,k3,k4;
    for(int i=1;i<n;i++){
        t[i]=i*step;

        x[0]=ic[0];
        v[0]=ic[1];

        k1=step*(-x[i-1]);
        k2=step*(-x[i-1]+step/2);
        k3=step*(-x[i-1]+step/2);
        k4=step*(-x[i-1]+step);
        v[i]=v[i-1]+(k1+2*k2+2*k3+k4)*step/6;

        k1=step*(v[i-1]);
        k2=step*(v[i-1]+step/2);
        k3=step*(v[i-1]+step/2);
        k4=step*(v[i-1]+step);
        x[i]=x[i-1]+(k1+2*k2+2*k3+k4)*step/6;
}}

void print_file(double* t, double* x, double* v, int n){
    int i;
    for(i=0; i<=n; i++)
    output<<setw(15)<<t[i]<<setw(15)<<x[i]<<setw(15)<<v[i]<<endl;
    return;}

Where are my errors?

Improve this question
$\endgroup$
2
$\begingroup$ 
    k1=step*(-x[i-1]);
    k2=step*(-x[i-1]+step/2);
    k3=step*(-x[i-1]+step/2);
    k4=step*(-x[i-1]+step);

That's not the form of RK4. Each k has a derivative calculation that uses the previous derivative estimate. The first should be step*f(x[i-1]), then the next should be step*f(x[i-1] + step/2 * k1), ... (written non-homogeneous). Then I don't see your derivative f defined anywhere, and then you also need to take that an make it work for a vector equation (probably make it inplace, etc.) since you have to update the x and v simultaneously. I would highly recommend just using software which is already made to do this since in the end, RK4 won't even be the most efficient way to handle this, even after you put in the work to somewhat test it.

Improve this answer
$\endgroup$
3
  • $\begingroup$ But in the case of the harmonic oscillator the f should be x, or not? This point confuses me. $\endgroup$ – mattiav27 Oct 29 '17 at 20:27
  • $\begingroup$ Well the system has two components: the position and the velocity. The derivative of the position is the velocity. The derivative of the velocity, for this model, is -x. So okay, I see you inlined the function (making it less general), but you didn't didn't add the k1 part to k2, etc. $\endgroup$ – Chris Rackauckas Oct 30 '17 at 2:20
  • $\begingroup$ Rackauakas Could you please edit your question with the correct version of the code? I will never learn in this way. $\endgroup$ – mattiav27 Oct 30 '17 at 7:53

Not the answer you're looking for? Browse other questions tagged or ask your own question.