Ошибка в больших итерациях n-мерного решателя Рунге-Кутты 4-го порядка

Код, представленный здесь, работает нормально, но, скажем, чтобы найти точки бифуркации при изменении омега для (-30,30), а не (10,30), таким образом, изменяя «int o» с 2000 на 6000, на экране появляется следующее сообщение:

Необработанное исключение в 0x7665B802 в Bifurcation_Plotter.exe: исключение Microsoft 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;
}