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