I am trying to solve a set of two 2nd order differential equations using the 4th order Runge-Kutta method. For this purpose, I have rewritten the two 2nd order equations into four 1st order equations, so that the problem gets reduced to solving four 1st order equations. I am writing the code in C language. The two variables in the problem are the radial coordinate (r) and the zenith angle (theta). I expect to obtain a curve like the following:
One can observe from the above plot that the derivative of theta changes at the turning points, but the derivative of r never changes its sign because the particle is always moving inwards towards the origin. So, the constraint is: the derivative of r is always positive.
I have written the following code where I find that soon after the derivative of theta changes sign, the derivative of r also changes sign and the particle starts to move in the outside direction, and I couldn't implement the constraint in the code. The plot that I obtained is:
The code that I have written is the following:
#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 = fabs(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=500.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,test;
FILE *fp1=NULL;
fp1=fopen("rthetaplot.txt","w");
// Integration
for(t=499.0;t>=0.0;t-=0.5)
{
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);
t0 = t;
r0 = r;
theta0 = theta;
R0 = R;
T0 = T;
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);
}
In the code, I have defined the 1st derivative of r and theta as R and T respectively. So, I have four 1st order differential equations for r, theta, R and T. I started to integrate from a point (r,theta,R,T)=(50,pi/3,0,0) and used a negative step-size.
So, can anyone please suggest how to incorporate the constraint that the derivative of r is always positive so that the particle continues to move inwards after crossing the turning points?
