c++常微分方程组

Question

您好:) 我尝试用四阶龙格-库塔方法求解 C++ 中常微分函数的方程组(一阶微分方程)。它们是两个耦合的 FitzHugh-Nagumo 模型神经元的 ODE，它描述了跨神经元细胞膜的静电势。通过运行此代码，我遇到了段错误。大家能帮帮我吗，我应该改什么代码？

float f1(float x1, float y1, float gamma1, float x2){
  float xx1;
  xx1 = pow(x1,3);
  return x1-(xx1/3)-y1+gamma1*x2;
}

float f2(float eps1, float x1, float a1) {
  return eps1 * (x1 +a1);
}

float f3(float x2, float y2, float gamma2, float x1) {
  return x2-(x2/3)-y2+gamma2*x1;
}

float f4(float eps2, float x2, float a2) {
  return eps2 * (x2 +a2);
}


int main(){

 // Set the parameters for the FitzHugh-Nagumo model
 // Either read them from the parameter.dat 
 // or simply hard code them here.
float eps1 = 0.1;
float gamma1 = 2.0;
float eps2 = 0.1;
float gamma2 = -1.5;

float a1[3] = {1.3, 1.3, 1.3};
float a2[3] = {1.3, 0.5, 0.25};

float LO = -2;
float HI = 2;

vector<float> fx1;
vector<float> fx2;
vector<float> fy1;
vector<float> fy2;

fx1[0] = LO + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(HI-LO)));
fx2[0] = LO + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(HI-LO)));

fy1[0] = LO + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(HI-LO)));
fy2[0] = LO + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(HI-LO)));

const int number_steps = 5000;
const float step_size = 0.05;

for(int i=0; i <= number_steps; ++i){

  float x1_1 = f1(fx1[i],fy1[i],gamma1,fx2[i]); 
  float y1_1 = f2(eps1, fx1[i], a1[0]); 
  float x2_1 = f3(fx2[i],fy2[i],gamma2,fx1[i]); 
  float y2_1 = f4(eps2,fx2[i],a2[0]); 

  float x1_2 = f1(fx1[i]+((step_size/2)*x1_1),fy1[i]+((step_size/2)*y1_1),gamma1,fx2[i]+((step_size/2)*x2_1)); 
  float y1_2 = f2(eps1, fx1[i]+((step_size/2)*x1_1), a1[0]); 
  float x2_2 = f3(fx2[i]+((step_size/2)*x2_1),fy2[i]+((step_size/2)*y2_1),gamma2,fx1[i]+((step_size/2)*x1_1)); 
  float y2_2 = f4(eps2,fx2[i]+((step_size/2)*x2_1),a2[0]); 

  float x1_3 = f1(fx1[i]+((step_size/2)*x1_2),fy1[i]+((step_size/2)*y1_2),gamma1,fx2[i]+((step_size/2)*x2_2)); 
  float y1_3 = f2(eps1, fx1[i]+((step_size/2)*x1_2), a1[0]); 
  float x2_3 = f3(fx2[i]+((step_size/2)*x2_2),fy2[i]+((step_size/2)*y2_2),gamma2,fx1[i]+((step_size/2)*x1_2)); 
  float y2_3 = f4(eps2,fx2[i]+((step_size/2)*x2_2),a2[0]); 

  float x1_4 = f1(fx1[i]+(step_size*x1_3),fy1[i]+(step_size*y1_3),gamma1,fx2[i]+(step_size*x2_3)); 
  float y1_4 = f2(eps1, fx1[i]+(step_size*x1_3), a1[0]); 
  float x2_4 = f3(fx2[i]+(step_size*x2_3),fy2[i]+(step_size*y2_3),gamma2,fx1[i]+(step_size*x1_3)); 
  float y2_4 = f4(eps2,fx2[i]+(step_size*x2_3),a2[0]); 

  fx1[i+1] = fx1[i] +(step_size/6)*(x1_1 + 2*x1_2 + 2*x1_3 + x1_4);
  fy1[i+1] = fy1[i] +(step_size/6)*(y1_1 + 2*y1_2 + 2*y1_3 + y1_4);
  fx2[i+1] = fx2[i] +(step_size/6)*(x2_1 + 2*x2_2 + 2*x2_3 + x2_4);
  fy2[i+1] = fy2[i] +(step_size/6)*(y2_1 + 2*y2_2 + 2*y2_3 + y2_4);

}

Answer 1

分配您的数据。

vector<float> fx1;
vector<float> fx2;
vector<float> fy1;
vector<float> fy2;

考虑到数据已经可用，创建您随后尝试访问的空 vector 。它不是。

使用:

vector<float> fx1(number_steps+2);
vector<float> fx2(number_steps+2);
vector<float> fy1(number_steps+2);
vector<float> fy2(number_steps+2);