#include <stdio.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <ctime>
using namespace std;
int tau_1=5, tau_2=50, g_slow=2;
double k_tau=0.05;
int g_fast=2;
double tau_m=0.16, I_app=0.697;
double step(double h,int n,double y[2], double r1);
double f(int i,double y[],double r1){
switch (i){
case 0:
return ((-y[0]+tanh(g_fast*y[0])-y[1]-I_app +0.1*r1)/tau_m);break;
case 1: return ((g_slow*y[0]-y[1])/(tau_2+(tau_1-tau_2)/(1+1/exp(y[0]/k_tau))));break;
default : break;
}
}
double step(double h,int n,double y[],double r1)
{
int i;
double yt[2];
double k1[2];
double k2[2];
double k3[2];
double k4[2];
for(i = 0; i <= n; i++)
{
y[i];
k1[i] = h*f(i,y,r1);
}
for(i = 0; i <= n; i++)
{
yt[i] = y[i]+0.5*k1[i];
}
for(i = 0; i <= n; i++)
{
k2[i] = h*f(i, yt,r1);
}
for(i = 0; i <= n; i++)
{
yt[i] = y[i]+0.5*k2[i];
}
for(i = 0; i <= n; i++)
{
k3[i] = h*f(i, yt,r1);
}
for(i = 0; i <= n; i++)
{
yt[i] = y[i]+k3[i];
}
for(i = 0; i <= n; i++)
{
k4[i] = h*f(i, yt,r1);
}
for(i = 0; i <= n; i++)
{
y[i] = y[i]+(k1[i]+2.0*k2[i]+2.0*k3[i]+k4[i])/6.0;
}
return y[2];
}
void spaik(double Y1[], double Y2[],double T[],int steps)
{
double ArrayTime[steps-1];
int j=0;
ofstream out2;
out2.open("Runge-Kutta2.txt");
for (int i=1;i<=steps-1;i++)
{
if ( Y1[i]>0.75 && Y1[i+1]>0.75 && Y1[i-1]>0.75)
{
if(Y1[i]>Y1[i+1] )
{
if(Y1[i]>Y1[i-1])
{
if(T[i])
{
ArrayTime[j]=T[i];
j=j+1;
out2<<T[i]<<"\t"<<endl;
}
}
}
}
}
double tauspike;
double ArrayTauSpike[steps-1];
for (int i1=0;i1<=j-2;i1++)
{
tauspike=ArrayTime[i1+1]-ArrayTime[i1];
ArrayTauSpike[i1]=tauspike;
cout<<tauspike<<"\t"<<endl;
}
}
void solvesystemrungekutta(double t,double t1,int steps,double y[]){
ofstream out;
out.open("Runge-Kutta.txt");
srand(time(0));
double r[steps-1];
double T[steps-1];
double Y1[steps-1];
double Y2[steps-1];
double r1;
int j=0;
for(int i = 1; i <= steps-1; i++)
{
y[2];
Y1[j]=y[0];
Y2[j]=y[1];
r[j]=1.0 * rand() / RAND_MAX;
r1=r[j];
step(((t1-t)/steps),2,y,r1);
T[j+1]=T[j]+((t1-t)/steps);
out<<T[j]<<"\t"<<Y1[j]<<"\t"<<Y2[j]<<endl;
j=j+1;
}
spaik(Y1,Y2,T,steps);
out.close();
}
int main(){
//первая координата
double temporaryresult[2]={0.1,0.1};
solvesystemrungekutta(1,500,10000,temporaryresult);
return 0;
}