Runge-Kutta 4th-Order

Question

Right now I am trying to make my Runge-Kutta work. Actually I have to solve a ODE-System which is a lot more complicated than this ODE-System, but the runge-kutta-function is the same. I'm solving this simple system to check if the simulation converges to the analytical solution.

First of all I would like to ask you if you see any mistakes in my code. The Problem is, that I get a solution but it differs from the analytical solution (solved it via symbolab). Here is my code:

#include <iostream>
#include <vector>
#include <fstream>

using namespace std;

double fx(double x, double y){
    return y;
}

double fy(double x, double y){
    return x + y;
}

void runge_kutta(double dt, double &x, double &y){
    vector<double> k(4), l(4);

    #define dt_half dt/2.0
    #define dt_six dt/6.0

    k.at(0)=fx(x, y);
    l.at(0)=fy(x,y);

    k.at(1)=fx(x+dt_half*k.at(0), y+dt_half*l.at(0));
    l.at(1)=fx(x+dt_half*l.at(0), y+dt_half*l.at(0));

    k.at(2)=fx(x+dt_half*k.at(1), y+dt_half*l.at(1));
    l.at(2)=fx(x+dt_half*l.at(1), y+dt_half*l.at(1));

    k.at(3)=fx(x+dt*k.at(2), y+dt*l.at(2));
    l.at(3)=fy(x+dt*k.at(2), y+dt*l.at(2));

    x+=dt_six*(k.at(0)+2.0*k.at(1)+2.0*k.at(2)+k.at(3));
    y+=dt_six*(l.at(0)+2.0*l.at(1)+2.0*l.at(2)+l.at(3));


}

double t=10.0;
double dt=1e-3;

double x=0.0;
double y=1.0;

ofstream out;

int main(){
    out.open("out.txt");
    out << "#t, x, y" << endl;
    for(double time=0.0;time<t;time+=dt){
        cout<<time<<endl;
        out << time << ' '
            << x << ' '
            << y << ' '
            << endl;
        runge_kutta(dt, x, y);
    }
    out.close();  
}

Can you see any errors?

I'm thankful in advance for any help!

Please add sample input, expected output, and real output from the program. — R Sahu, Aug 20 '21 at 19:58
Your `l_2` and `l_3` are using the wrong function, are they not? Also, why do you even have an `l` in the first place? — Ext3h, Aug 20 '21 at 20:03
_"Debug my code for me"_ is off-topic here. Please read [How to debug small programs.](//ericlippert.com/2014/03/05/how-to-debug-small-programs/) and use a [debugger](//stackoverflow.com/q/25385173/843953). You need to perform basic diagnosis to include with your post. At the very least, print the suspected values at the point of error and trace them back to their sources Show where the intermediate results deviate from the ones you expect. — Pranav Hosangadi, Aug 20 '21 at 21:35

score 1 · Answer 1 · answered Aug 20 '21 at 21:32

1

Look closely at these lines

l.at(1)=fx(x+dt_half*l.at(0), y+dt_half*l.at(0));
l.at(2)=fx(x+dt_half*l.at(1), y+dt_half*l.at(1));

they should use fy instead of fx I think.

answered Aug 20 '21 at 21:32

Pepijn Kramer

9,356
2
8
19

score 1 · Accepted Answer · answered Aug 21 '21 at 00:50

1

I have highlighted several inconsistencies in the RK4 code

Anything that deals with y should call fy and have derivative l and anything that deals with x has derivative k.

answered Aug 21 '21 at 00:50

John Alexiou

28,472
11
77
133

Runge-Kutta 4th-Order

2 Answers2