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?

$\endgroup$

closed as off-topic by Chris Rackauckas, Kirill, Christian Clason, Wolfgang Bangerth, nicoguaro Oct 30 '17 at 3:32

This question appears to be off-topic. The users who voted to close gave this specific reason:

If this question can be reworded to fit the rules in the help center, please edit the question.

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.

$\endgroup$
  • $\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.