如何在四阶Runge-Kutta法中添加r导数恒正的约束？

问题 #

我用四阶Runge-Kutta（RK4）方法求解一组二阶微分方程，已将其转化为四个一阶微分方程，用C语言实现。变量定义如下：

• 径向坐标r，天顶角theta
• R = dr/dt，T = dtheta/dt

根据物理模型，粒子应始终向内运动，即R需保持恒正，但当前代码中T变号后，R也随之变号，粒子开始向外运动，无法得到预期轨迹。初始积分点为(r,theta,R,T)=(50,π/3,0,0)，采用负步长（h=-0.001）积分，希望在RK4算法中加入R恒正的约束。

当前代码：

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

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

// dr/dt 的被积函数
double integrand1(double R) 
{ 
    return R;
}

// dtheta/dt 的被积函数
double integrand2(double T) 
{ 
    return T;
}

// 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;
}

// 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;
}


// -------------------- RUNGE-KUTTA 算法 ----------------------------

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;
} 

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

// ------------------
// 主函数
//-------------------

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");


    // 积分循环
    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);
} 

预期轨迹：粒子持续向内螺旋运动，轨迹收敛于中心区域
实际轨迹：粒子运动一段时间后开始向外发散，偏离预期路径

解决方案 #

1. 修正物理模型的符号一致性 #

首先确认积分方向与R符号的逻辑：
采用负步长h=-0.001时，时间t从500递减到0，粒子向内运动意味着径向坐标r随时间减小。由于R=dr/dt，dr = R*dt，而dt=h<0，要让dr<0（r减小），R必须为正，这个逻辑是对的。

但需要检查dR/dt的被积函数integrand3的符号是否符合物理规律：
term1 = r*T²是离心力项，会阻碍粒子向内运动，对应dR/dt应该为负，但当前代码中是+term1，这可能是导致R变号的根源。建议修正该项符号：

term1 = -(r*T*T); // 改为负号，符合离心力阻碍向内运动的物理规律

2. 在RK4迭代中加入R的强制约束 #

如果确认微分方程正确，仅需在数值积分中保证R恒正，可在RK4每一步更新R后加入约束：

// 在RK4循环内更新R后添加
R = R + (1.0/6.0)*(m1 + 2*m2 + 2*m3 + m4);
if (R < 1e-8) {
    R = 1e-8; // 限制R为极小正值，避免变为负或0导致积分停滞
}

也可在计算中间步（k2/k3/k4）时提前约束R的中间值，防止中间计算出现负号影响后续迭代：

// 比如计算k2时
double R_mid = R + 0.5*m1;
if (R_mid < 1e-8) R_mid = 1e-8;
k2 = h*integrand1(R_mid);

3. 调整积分步长或改用自适应步长 #

当前固定步长h=-0.001可能过大，导致数值振荡进而使R变号。可尝试减小步长（比如h=-0.0001），或改用自适应步长RK4算法，根据每一步的误差动态调整步长，提升数值稳定性。

修正后的RK4核心部分示例 #

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);

    double R_mid1 = R + 0.5*m1;
    if (R_mid1 < 1e-8) R_mid1 = 1e-8;
    k2 = h*integrand1(R_mid1);
    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_mid1, T + 0.5*n1);
    
    double R_mid2 = R + 0.5*m2;
    if (R_mid2 < 1e-8) R_mid2 = 1e-8;
    k3 = h*integrand1(R_mid2);
    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_mid2, T + 0.5*n2);

    double R_mid3 = R + m3;
    if (R_mid3 < 1e-8) R_mid3 = 1e-8;
    k4 = h*integrand1(R_mid3);
    l4 = h*integrand2(T + n3);
    m4 = h*integrand3(r + k3, theta + l3, T + n3);
    n4 = h*integrand4(r + k3, theta + l3, R_mid3, 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);
    if (R < 1e-8) R = 1e-8; // 最终约束
    T = T + (1.0/6.0)*(n1 + 2*n2 + 2*n3 + n4);      

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

内容的提问来源于stack exchange，提问作者AnnaD