n维四阶Runge-Kutta求解器的大迭代中的误差

Question

这里给出的代码工作正常，但是说当找到（-30,30）而不是（10,30）的ω时发现分叉点，从而改变'INT O' 从2000到6000在屏幕上显示以下消息，在在0x7665B802在Bifurcation_Plotter.exe

未处理的异常：微软C++异常：则为std :: bad_alloc在存储器位置0x012FF544。

时间步骤需要保持原样以确保结果的准确性。

所有帮助非常感谢:)

//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个双打，而不是动态分配一些结构，保护访问，使用复制操作，我认为你的问题是围绕动态分配。

我看到一些内循环被称为大约10万次（参见o（2K）和n（5K）的顺序）。函数r4kvec返回一个指向永不释放的动态分配区域的指针。电话一样）因此10MLN * 6 * 64，让我觉得你可能会很接近exaust内存（但也要看你正在运行它的系统）

所以说：

• delete[]每次拨打r4kvec返回的数据
• ，因为new通常比6双倍的副本更昂贵，所以考虑使用6双精心设计的结构以便使用堆栈并降低复杂性。

希望这可以帮助， 斯特凡诺

---

所以像这样的循环：

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