Numerically solve a system of ODEs using 4th order Runge-Kutta method integrating backwards

numerical methodsordinary differential equationsrunge-kutta-methodssystems of equations

I am trying to solve a system of two 2nd order ODEs using the 4th order Runge-Kutta (RK4) method. The equations are of the form:
$$\frac{d^2r}{dt^2}=f(r,\theta,\dot{r},\dot{\theta}),$$
$$\frac{d^2\theta}{dt^2}=g(r,\theta,\dot{r},\dot{\theta}),$$
where the overdot means time-derivative, i.e., $\dot{x}=\dfrac{dx}{dt},$ for any variable $x$.

Now, we can rewrite the above equations in the following form:
$$\frac{dr}{dt}=R,$$
$$\frac{d\theta}{dt}=\Theta,$$
$$\frac{dR}{dt}=f(r,\theta,R,\Theta),$$
$$\frac{d\Theta}{dt}=g(r,\theta,R,\Theta).$$

These 4 equations can be easily solved using the RK4 algorithm. Also, note that, none of the integrands depend explicitly on the variable $t$.

Problem:

I need to solve these four 1st order ODEs and plot the $r-\theta$ trajectory in the backward direction. This means that the initial values of $(r,\theta)$ are provided for a large value of $r$ and then the integration is performed backwards towards the origin. But, I am having the problem regarding how to implement this in the RK4 algorithm. What changes needs to be done in the RK4 code written for the forward integration?

I am attaching the RK4 code as follows:

#include<stdio.h>
#include <math.h>

// ***********************************************************************************************

// The integrand for dr/dt 

double integrand1(double R) 
{ 
    return R;
}

// The integrand for dtheta/dt 

double integrand2(double T) 
{ 
    return T;
}

// The integrand for dR/dt 

double integrand3(double r, double theta, double T)
{
    double rnot, a=0.8, l=2.68, r1=3.5, term1, term2, term3, term4, term5, integrand3;       
    rnot = pow(a,0.125);
    term1 = (r*T*T);
    term2 = (1/((r-rnot)*(r-rnot)));
    term3 = ((l*l)/(2*r*r*r*sin(theta)*sin(theta)));
    term4 = ((3*(r-2)*l*l)/(2*r*r*r*r*sin(theta)*sin(theta)));
    term5 = ((3*r1*a*l)/(r*r*r*r));
    integrand3 = term1 - term2 - term3 + term4 + term5;
    return integrand3;
}

// The integrand for dT/dt 

double integrand4(double r, double theta, double R, double T)
{
    double l=2.68, term1, term2, integrand4;                                               
    term1 = ((r-2)*l*l*cos(theta))/(r*r*r*r*r*sin(theta)*sin(theta)*sin(theta));
    term2 = ((2*R*T)/r);
    integrand4 = term1 - term2;
    return integrand4;
}


// --------------- THE RUNGE-KUTTA ALGORITHM --------------------

double rk4(double t0, double r0, double theta0, double R0, double T0, double t, double h, double arr[])    
{ 
    int n = (int)(fabs((t-t0)/h)); 
    double k1,k2,k3,k4,l1,l2,l3,l4,m1,m2,m3,m4,n1,n2,n3,n4;

    double r = r0;
    double theta = theta0;
    double R = R0;
    double T = T0;

    for (int i=1; i<=n; i++) 
        {
        k1 = h*integrand1(R);
        l1 = h*integrand2(T);
        m1 = h*integrand3(r, theta, T);
        n1 = h*integrand4(r, theta, R, T);

        k2 = h*integrand1(R + 0.5*m1);
        l2 = h*integrand2(T + 0.5*n1);
        m2 = h*integrand3(r + 0.5*k1, theta + 0.5*l1, T + 0.5*n1);
        n2 = h*integrand4(r + 0.5*k1, theta + 0.5*l1, R + 0.5*m1, T + 0.5*n1);
        
        k3 = h*integrand1(R + 0.5*m2);
        l3 = h*integrand2(T + 0.5*n2);
        m3 = h*integrand3(r + 0.5*k2, theta + 0.5*l2, T + 0.5*n2);
        n3 = h*integrand4(r + 0.5*k2, theta + 0.5*l2, R + 0.5*m2, T + 0.5*n2);

        k4 = h*integrand1(R + m3);
        l4 = h*integrand2(T + n3);
        m4 = h*integrand3(r + k3, theta + l3, T + n3);
        n4 = h*integrand4(r + k3, theta + l3, R + m3, T + n3);

        r = r + (1.0/6.0)*(k1 + 2*k2 + 2*k3 + k4);
        theta = theta + (1.0/6.0)*(l1 + 2*l2 + 2*l3 + l4);
        R = R + (1.0/6.0)*(m1 + 2*m2 + 2*m3 + m4);
        T = T + (1.0/6.0)*(n1 + 2*n2 + 2*n3 + n4);      

        arr[0] = r;     
        arr[1] = theta;
        arr[2] = R;
        arr[3] = T;

        t0 = t0 + h; 
    } return 0;
} 



// ------------------------------------------------------------------------------------

// ------------------
// The main function
//-------------------

int main() 
{ 
double t0=1000.0,r0=50.0,theta0=1.0467,R0=0.0,T0=0.0,t,r,theta,R,T,a,h=-0.001,x,z,arr[4],i;

FILE *fp1=NULL;
fp1=fopen("rthetaplot.txt","w");


// Integration

for(t=999.0;t>=0.0;t-=1.0)
{
rk4(t0,r0,theta0,R0,T0,t,h,arr);
r = arr[0];
theta = arr[1];
R = arr[2]; 
T = arr[3];
x = r*sin(theta);
z = r*cos(theta);
fprintf(fp1,"%f \t %f \t %f \t %f \t %f \t %f \t %f \n", t,r,theta,R,T,x,z);
}


fclose(fp1);
}

Can anyone please suggest what necessary changes are required so that the backward integration can be easily performed?

EDIT: I have entered the entire code.

Best Answer

You need to change nothing. Just stay with $(t_0,u_0)$ as the initial point where the initial condition is given, $t$ or $t_{end}$ the point where the integration has to end, and use a negative step size $h$ to move backwards in time if $t_{end}<t_0$.

Related Question