解決微分方程式は

Question

主機能は、メインループ反復に計算をフィードバックgate_probabilities別の関数に見出さ順次相互依存微分方程式のセリエにベクトル値を供給するループを実行します。

方程式の解は、対応するベクトルにそれぞれ連続して入力されます。これは、両方の関数で行われます。中間式はm_new; h_new; n_newの場合、時間に対して最終的にプロットされる主ベクトルはVです。 （コードではないプロット）

ベクトルが関数本体の外で定義されている（これらはグローバルであると考えるが、私は、問題は、彼らが適切にそのように宣言していないと思います）

ビルには、エラーを生成しません、私は実行時にデバッグエラーが発生します： "Abort（）が呼び出されました"。 timevectorが空であるよう

は、任意の助けをいただければ幸いです:)

#include "C:\Users\jenbe\Documents\Visual Studio 2015\std_lib_facilities.h" 
#include <math.h> 
#include <iostream> 
#include <cmath> 
#include <tuple> 
#include <vector> 

using namespace std; 

vector<double> mvector; 
vector<double> hvector; 
vector<double> nvector; 
vector<double> V; // potential vector 
vector<double> timevector; 

//global constants 

double gL = 0.1; // mS/cm^2 
double gK = 9; 
double gNa = 35; 
double EL = -65; // mV 
double EK = -90; 
double ENa = 55; 
double phi = 5; // Coefficient increasing reaction speed of the alpha and beta constants which it multiplies. 
double Iapp = 5; // microAmpere values (as in Wang-Buzsaki) 
double runtime = 3000;// 3 Seconds 
double dt = 0.001; // Timestep = 1ms 
double V_init = -60; 


std::tuple<double> gate_probabilities(double v, double m, double h, double n) { 
    double am, bm, ah, bh, an, bn, dh, dn; 

    //first stage -> gate probabilities 

    am = -0.1*(v + 35)/(exp(-0.1*(v + 35)) - 1); //the probability of a closed gate to open 
    bm = 4 * exp(-(v + 60)/18); //the probability of an open gate to be closed. 
    ah = 0.07*exp(-(v + 58)/20); 
    bh = 1/(exp(-0.1*(v + 28)) + 1); 
    an = -0.01*(v + 34)/(exp(-0.1*(v + 34)) - 1); 
    bn = 0.125 * exp(-(v + 44)/80); 

    //second stage -> gate states 

    dh = phi * (ah*(1 - h) - bh*h); // inactivation variable h 
    dn = phi * (an*(1 - n) - bn*n); // inward recitifier 

    double m_new = am/(am + bm); //activation variable 
    double h_new = dh*dt + h; 
    double n_new = dn*dt + n; 

    mvector.push_back(m_new); 
    hvector.push_back(h_new); 
    nvector.push_back(n_new); 

    // third stage -> currents 

    double IL = gL*(v - EL); 
    double INa = gNa * pow(m, 3) * h * (v - ENa); 
    double IK = gK * pow(n, 4) * (v - EK); // delayed recitifier 
    double currents = -INa - IK - IL + Iapp; 

    return std::make_tuple(currents); 
} 


int main() { 

      //initialize vectors - later modify with function within function 
    double ami = -0.1*(V_init + 35)/(exp(-0.1*(V_init + 35)) - 1); 
    double bmi = 4 * exp(-(V_init + 60)/18); 
    double ahi = 0.07*exp(-(V_init + 58)/20); // 
    double bhi = 1/(exp(-0.1*(V_init + 28)) + 1); 
    double ani = -0.01*(V_init + 34)/(exp(-0.1*(V_init + 34)) - 1); 
    double bni = 0.125 * exp(-(V_init + 44)/80); 
    V.push_back(V_init); //V_init 
    double m_init = (ami/(ami + bmi)); 
    double h_init = (ahi/(ahi + bhi)); 
    double n_init = (ani/(ani + bni)); 

    mvector.push_back(m_init); 
    hvector.push_back(h_init); 
    nvector.push_back(n_init); 

    for (int i{ 1 }; i <= runtime; ++i) { // 1ms iterations over 3000ms range 

     auto ret = gate_probabilities(V[i - 1], mvector[i - 1], hvector[i - 1], nvector[i - 1]); 
     double currents; 
     currents = std::get<0>(ret); 

     double potential = currents*dt + V[i - 1]; //calculate new V. 
     V.push_back(potential); //incorporates new V into a vector. 
     timevector.push_back(timevector[i - 1] + dt); 
    } 
} 

Answer 1

をあなたのライン

timevector.push_back(timevector[i - 1] + dt); 

は、非常に最初のループの実行では動作しません。