n维四阶龙格-库塔求解器大迭代中的误差

Question

此处提供的代码可以正常工作，但要在 (-30,30) 而不是 (10,30) 改变 omega 时找到分叉点，因此将“int o”从 2000 更改为 6000，屏幕上会出现以下消息,

Bifurcation_Plotter.exe 中 0x7665B802 处的未处理异常：Microsoft C++ 异常：内存位置 0x012FF544 处的 std::bad_alloc。

时间步长需要保持不变以确保结果的准确性。

非常感谢所有帮助:)

//NOTE: this code has memory issues, if compiling be careful to adjust step size to obtain the desired plot
    //This program computes solutions for I_A or I_B and stores them to an array
    //This array is then evaluated to find the local maxima
    //Further evaluation finds what the local maxima settle to
    //These settled values are found for a varying omega to produce a bifurcation plot


    #include <cstdlib>
    #include <iostream>
    #include <string>
    #include <fstream>
    #include <iomanip>
    #include <cmath>

        using namespace std;

        double *CoupledLaser_rhs(double t, int m, double x[]);
        double *rk4vec (double t0, int m, double u0[], double dt, double *f(double t, int m, double x[]));

        //Global parameter values
            double beta = 8.5;
            double gama = 10.0;
            double alpha = 2.0;
            double kappa = 39.97501809;//d=1.2
            double omega;
            double lambda = 2;

        int main()
        {
            //file to store bif points
            string bif_points_filename = "(10,30)mod_bif_points_IA_d=1.2.txt";
            ofstream bif_unit;


            //


            double dt = 0.01;//time-step
            double domega = 0.01;//omega-step
            int i;
            int j;
            int k;
            int l;
            int m = 6;//no. of dimensions
            int n = 5000;//no. of time evaluation steps (n*dt = time)
            int o = 2000;//no. of delta_omega evaluation steps
            double t;
            double *x;
            double *xnew;
            double current;

            //

            cout<<"\n";
            cout<<"CoupledLaser_RKSolver\n";
            cout<<"Compute solutions of the Coupled Laser system.\n";
            cout<<"Write data to file.\n";

            //


            //I.C.'s in 0th entry//
            t = 0.0;
            omega=10;

            x = new double [m];
            x[0] = 1.0;
            x[1] = 1.0;
            x[2] = 1.0;
            x[3] = 1.0;
            x[4] = 0.001;
            x[5] = 0.001;


            //


            //define array to store elements of I_A
            double *arr = new double[1000];

        //


        //Approximate solution at equally spaced times of time step dt//
        bif_unit.open(bif_points_filename.c_str());
        for(l=0; l<o; l++)
        {
            for(j=0; j<n-1000; j++)
                {
                    current = ((x[0])*(x[0])+(x[1])*(x[1]));
                    xnew = rk4vec(t, m, x, dt, CoupledLaser_rhs);
                    for(i=0; i<m; i++)
                    {
                        x[i] = xnew[i];
                    }
                    t=t+dt;
                }   
            arr[0]=current;
            for(j=0; j<1000; j++)
                {
                    arr[j]=((x[0])*(x[0])+(x[1])*(x[1]));
                    xnew = rk4vec(t, m, x, dt, CoupledLaser_rhs);
                    for(i=0; i<m; i++)
                    {
                        x[i] = xnew[i];
                    }
                    t=t+dt;
                }
            for(k=50; k<1000-50; k++)
                {
                    if(arr[k]>arr[k+1] && arr[k]>arr[k-1])
                        {
                            bif_unit <<omega<<","<< arr[k]<<"\n";
                        }
                }
            omega = omega + domega;
        }

        bif_unit.close();

        //

    cout << "Created local maxima vs omega bifurcation file " << bif_points_filename<<"\".\n";



        //END//

        cout<<"\n";
        cout<<"CoupledLaser_ODE:\n";
        cout<<"Normal end of execution.\n";
        cout<<"\n";
    }

//

//Evaluates the rhs of the coupled laser field equations
//t; value of the independent time variable, m; spatial dimension, x[]; values of the dependent variables at time t
//Output; values of the derivatives of the dependent variables at time t
//x[0] = E_Ax, x[1] = E_Ay, x[2] = E_Bx, x[3] = E_By, x[4] = N_A, x[5] = N_B

double *CoupledLaser_rhs(double t, int m, double x[])
{
    double *dxdt;
    dxdt = new double [m];

    dxdt[0] = beta*gama*(x[4]*x[0]) + alpha*beta*gama*(x[4]*x[1]) - kappa*x[3];
    dxdt[1] = beta*gama*(x[4]*x[1]) - alpha*beta*gama*(x[4]*x[0]) + kappa*x[2];
    dxdt[2] = beta*gama*(x[5]*x[2]) + alpha*beta*gama*(x[5]*x[3]) - kappa*x[1] + omega*x[3];
    dxdt[3] = beta*gama*(x[5]*x[3]) - alpha*beta*gama*(x[5]*x[2]) + kappa*x[0] - omega*x[2];
    dxdt[4] = lambda - x[4] - 1 - (x[0])*(x[0]) - (x[1])*(x[1]) - beta*(x[4])*((x[0])*(x[0])+(x[1])*(x[1]));
    dxdt[5] = lambda - x[5] - 1 - (x[2])*(x[2]) - (x[3])*(x[3]) - beta*(x[5])*((x[2])*(x[2])+(x[3])*(x[3]));



    return dxdt;
}

//IVP of the form du/dt = f(t,u) & u(t0) = u0
//User supplies the current values of t, u, step-size dt, and a function to evaluate the derivative, the function can compute the 4th-order Runge-Kutta estimate to the solution at time t+dt
//t0; current time, m; dimension of space, u0[]; solution estimate at current time, dt: time-step, *f; function which evaluates the derivative of the rhs of problem
//Output; 4th-order Runge-Kutta solution estimate at time t0+dt

double *rk4vec(double t0, int m, double x0[], double dt, double *f(double t, int m, double x[]))
{
    double *k1;
    double *k2;
    double *k3;
    double *k4;
    double t;
    double *x1;
    double *x2;
    double *x3;
    int i;
    double *x;


    //four sample values of the derivative

    k1 = f(t0, m, x0);

    t = t0 + dt/2.0;
    x1 = new double[m];

    for(i=0; i<m; i++)
    {
        x1[i] = x0[i] + dt*(k1[i]/2.0);
    }
    k2 = f(t, m, x1);

    x2 = new double[m];

    for(i=0; i<m; i++)
    {
        x2[i] = x0[i] + dt*(k2[i]/2.0);
    }
    k3 = f(t, m, x2);

    x3 = new double[m];
    for(i=0; i<m; i++)
    {
        x3[i] = x0[i] + dt*k3[i];
    }
    k4 = f(t0 + dt, m, x3);

    //combine to estimate solution

    x = new double[m];
    for(i=0; i<m; i++)
    {
        x[i] = x0[i] + dt*(k1[i] + 2.0*(k2[i]) + 2.0*(k3[i]) + k4[i])/(6.0);
    }

    //free memory

    delete [] k1;
    delete [] k2;
    delete [] k3;
    delete [] k4;
    delete [] x1;
    delete [] x2;
    delete [] x3;

    return x;
}

Answer 1

乍一看，因为我认为代码可以更好地编写，使用一些 6 个双精度的结构而不是动态分配，保护访问，使用复制运算符，我认为你的问题是动态分配。

我看到一些内部循环被调用了大约 1000 万次（参见 o (2K) 和 n (5K) 的顺序。函数r4kvec 返回一个指向永远不会释放的动态分配区域的指针（没有调用确实如此）。因此 10MLN * 6 * 64bit 让我认为您可能非常接近耗尽内存（但这也取决于您运行它的系统）。

这么说：

• delete[] 每个 r4kvec 调用返回的数据。
• 由于new 通常比 6 double 的副本更昂贵，因此请考虑设计一个精心设计的 6 double 结构以使用堆栈并降低复杂性。

希望这可能会有所帮助， 斯特凡诺

---

这样的循环：

    for(j=0; j<n-1000; j++)
        {
            current = ((x[0])*(x[0])+(x[1])*(x[1]));
            xnew = rk4vec(t, m, x, dt, CoupledLaser_rhs);
            for(i=0; i<m; i++)
            {
                x[i] = xnew[i];
            }
            t=t+dt;
        }  

应该变成：

    for(j=0; j<n-1000; j++)
        {
            current = ((x[0])*(x[0])+(x[1])*(x[1]));
            xnew = rk4vec(t, m, x, dt, CoupledLaser_rhs);
            for(i=0; i<m; i++)
            {
                x[i] = xnew[i];
            }
            t=t+dt;
            delete[] xnew; //release it!
        }