C ++ Система обыкновенных дифференциальных уравнений

A.SC спросил: 03 ноября 2018 в 08:43 в: c++

Привет :) Я пытаюсь решить систему уравнений обыкновенных дифференциальных функций в C ++ (дифференциальные уравнения первого порядка) методом Рунге-Кутты четвертого порядка. Они являются ODE двух связанных нейронов модели FitzHugh-Nagumo, они описывают электростатический потенциал через клеточную мембрану нейрона. При запуске этого кода я получаю ошибку сегментации. Может ли каждый помочь мне, что я должен изменить в коде?

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);}

0 ответов